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ABSTRACT 

We present a general formulation of special-relativistic magnetohydrodynamics and derive exact radially self- 
similar solutions for axisymmetric outflows from strongly magnetized, rotating compact objects. We generalize 
previous work by including thermal effects and analyze in detail the various forces that guide, accelerate, and 
collimate the flow. We demonstrate that, under the assumptions of a quasi-steady poloidal magnetic field and of 
a highly relativistic poloidal velocity, the equations become effectively time-independent and the motion can be 
described as a frozen pulse. We concentrate on trans-Alfvenic solutions and consider outflows that are super- 
Alfvenic throughout in the companion paper. Our results are applicable to relativistic jets in gamma-ray burst 
(GRB) sources, active galactic nuclei, and microquasars, but our discussion focuses on GRBs. We envision the 
outflows in this case to initially consist of a hot and optically thick mixture of baryons, electron-positron pairs, and 
photons. We show that the flow is at first accelerated thermally but that the bulk of the acceleration is magnetic, 
with the asymptotic Lorentz factor corresponding to a rough equipartition between the Poynting and kinetic-energy 
fluxes (i.e., ~ 50% of the injected total energy is converted into baryonic kinetic energy). The electromagnetic 
forces also strongly collimate the flow, giving rise to an asymptotically cylindrical structure. 

Subject headings: galaxies: jets — gamma rays: bursts — ISM: jets and outflows — MHD — methods: 
analytical — relativity 



1. INTRODUCTION 

The acceleration and collimation of powerful bipolar out- 
flows and jets in a variety of astronomical settings are often 
attributed to the action of magnetic fields (see, e.g., Livio 2000 
and Konigl & Pudritz 2000 for reviews). The commonly in- 
voked scenario is that magnetic field lines threading a rotating 
compact object or its surrounding accretion disk can efficiently 
tap the rotational energy of the source and accelerate gas to su- 
personic speeds through centrifugal and/or magnetic pressure- 
gradient forces. It is argued that the hoop stresses of the twisted 
field lines can account for the narrowness of many jets and 
that, in many cases, alternative production mechanisms (such 
as thermal driving) can be excluded on observational grounds. 

Although numerical simulations have provided useful in- 
sights into various aspects of hydromagnetic jet production, 
practical limitations have necessitated complementing this ap- 
proach with analytic studies. Owing to the complexity of the 
problem, the most general semianalytic solutions obtained so 
far have been time-independent and self-similar, patterned af- 
ter the pioneering disk-outflow solutions of Blandford & Payne 
(1982). The advantage of pursuing such solutions is that they 
are exact and self-consistent, and that they can be systemat- 
ically classified (e.g., Vlahakis & Tsinganos 1998). Further- 
more, these solutions are evidently rich enough to capture most 
of the relevant physics, as corroborated by numerical calcu- 
lations (e.g., Ouyed & Pudritz 1997; Ustyugova et al. 1999; 
Krasnopolsky, Li, & Blandford 1999). 

Almost all of the previous semianalytic work on jet magne- 
tohydrodynamics (MHD) was done in the Newtonian limit of 
nonrelativistic bulk and random speeds. Relativistic outflows 
are, however, observed quite commonly in Nature — with ex- 
amples including active galactic nuclei (AGNs), Galactic su- 
perluminal sources (often referred to as "microquasars"), pul- 



sars, and gamma-ray bursts (GRBs) — and in many of these 
cases magnetic fields are again implicated as the main driving 
and collimation mechanism (e.g., Blandford 2002a). This pro- 
vides a motivation for generalizing the Newtonian self-similar 
outflow solutions, although it is readily seen that this cannot be 
done in a totally straightforward manner. For one thing, special- 
relativistic MHD (unlike the nonrelativistic theory) involves a 
characteristic speed (the speed of light c), which precludes the 
incorporation of gravity into the self-similar equations and a 
simple matching of the outflow solution to a particular (e.g., 
Keplerian) disk rotation law. Furthermore, again in contrast 
with the nonrelativistic formulation, the displacement current 
and the charge density cannot be neglected in the constitutive 
equations (which now must also satisfy relativistic covariance). 

Despite the aforementioned complications, Li, Chiueh, & 
Begelman (1992) and Contopoulos (1994) succeeded in gen- 
eralizing the "cold," radially self-similar solutions of Bland- 
ford & Payne (1982) to the relativistic regime. Their solutions 
are characterized by the thermal pressure playing a negligi- 
ble role in the flow acceleration and by the flow being trans- 
Alfvenic: the poloidal speed is less than the poloidal compo- 
nent of the Alfven speed at the base of the flow, and comes to 
exceed it further up. Our aim in this paper is to further gen- 
eralize these solutions to the "hot" relativistic case — i.e., we 
allow both the bulk and the random speeds to be relativistic. 
This is motivated primarily by the desire to apply these solu- 
tions to GRB outflows, in which thermal driving by an opti- 
cally thick, hot "fireball" composed primarily of radiation and 
electron-positron pairs could play a role. In fact, most previ- 
ous models of GRB outflows were purely hydrodynamical and 
included only thermal driving (see, e.g., Piran 1999 for a re- 
view). It has subsequently been realized that energy deposition 
by annihilating neutrinos, which had been one of the main pro- 
posed heating mechanisms, would typically be inefficient (e.g., 
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Di Matteo, Perna, & Narayan 2002), and it was, furthermore, 
argued (Daigne & Mochkovitch 2002) that only a small frac- 
tion (< 1%) of the energy deposited in the source could be in 
thermal form to avoid generating strong photospheric emission 
in the outflow (for which there has been no observational ev- 
idence). For these (and other) reasons, it is now believed that 
magnetic fields play the dominant role in the driving of GRB 
jets (e.g., Meszaros 2002). Nevertheless, some thermal energy 
injection, either by neutrinos or by magnetic energy dissipation, 
is likely to take place (e.g., Narayan, Paczyhski, & Piran 1992; 
Meszaros, Laguna, & Rees 1993), and, as we show in §4.1, 
even if it contributes only a small fraction of the initial energy 
flux in the flow, it can dominate the early phase of the accel- 
eration. There are indications that thermal energy deposition 
at the source may also contribute to the initial acceleration of 
relativistic jets in AGNs (e.g., Melia, Liu, & Fatuzzo 2002; N. 
Vlahakis & A. Konigl, in preparation). It thus appears that, to 
fully understand the nature of such jets, it is necessary to model 
them in terms of a "hot" MHD outflow. 

The ability of large-scale, ordered magnetic fields to guide, 
collimate, and accelerate relativistic outflows has been previ- 
ously discussed by various authors (including, in the GRB con- 
text, Usov 1994, Thompson 1994; Meszaros & Rees 1997, 
Katz 1997, Kluzniak & Ruderman 1998, Lyutikov & Blackman 
2001, and Drenkhahn & Spruit 2002). These discussions, how- 
ever, did not include exact solutions from which detailed quan- 
titative estimates could be made. In a previous paper, Vlahakis 
& Konigl (200 1 , hereafter VK0 1 ) presented a semianaly tic self- 
similar solution of the "hot" relativistic MHD equations and 
applied it to the interpretation of GRBs. The full formalism un- 
derlying this solution is described in the present paper, where 
we also analyze its dependence on the relevant physical param- 
eters and compare it with other characteristic solutions of trans- 
Alfvenic flows. This discussion is extended in the companion 
paper (Vlahakis & Konigl 2003, hereafter Paper II), where we 
focus on super- Alfvenic flows. 

The self-similar solutions that we derive correspond to or- 
dered magnetic field configurations in the ideal-MHD limit. 
Although it is quite conceivable that the fields that drive the 
flow from a differentially rotating star or disk are at least in 
part small-scale and disordered, the statistical (temporal and 
spatial) averages of such fields could in principle have a sim- 
ilar effect to that of large-scale, ordered field configurations 
in providing both acceleration (Heinz & Begelman 2000) and 
collimation (Li 2002). The applicability of ideal MHD to the 
acceleration of ultrarelativistic flows has been questioned by 
Blandford (2002b; see also Lyutikov & Blandford 2002), who 
proposed instead a force-free electromagnetic formulation. We 
note in this connection that a force-free behavior can be recov- 
ered from the relativistic MHD formulation as a limiting case 
of negligible particle inertia. Furthermore, electromagnetic en- 
ergy dissipation — perhaps the strongest argument against the 
ideal-MHD modeling framework — has been claimed to lead, 
on its own, to an efficient conversion of Poynting flux into a 
highly relativistic bulk motion (e.g., Drenkhahn & Spruit 2002). 
Our ideal-MHD solutions may thus be regarded as a first step 
toward a more comprehensive theory in which dissipation ef- 
fects will be taken into account (and possibly further enhance 
the effectiveness of the magnetic acceleration process). 

Although the formulation presented in this paper is quite gen- 
eral, the main application that we consider is to GRB sources. 
For definiteness, we adopt the "internal shock" scenario for the 



origin of the prompt high-energy emission in GRBs (e.g., Pi- 
ran 1999). Considerations involving variability time scales (as 
interpreted in the context of this picture) as well as source en- 
ergetics support the identification of an accretion disk around 
a newly formed black hole as the source of the GRB outflow 
(e.g., Piran 2001a); accordingly, we concentrate on modeling 
jets from accretion disks. However, our solutions should be 
at least qualitatively applicable also to other source configu- 
rations, such as a rapidly spinning neutron star (e.g., Usov 
1994; Kluzniak & Ruderman 1998) or a rotating, magnetically 
threaded black hole (e.g., Blandford & Znajek 1977; van Put- 
ten 2001). Since GRB outflows have a limited duration (the 
value of which is plausibly related to the disk accretion time), 
a naive application of a steady-state similarity solution is not 
warranted. In previous, purely hydrodynamical models of GRB 
outflows, this difficulty was circumvented by applying the so- 
called "frozen pulse" approximation (Piran, Shemi, & Narayan 
1993). In this paper we prove (§2.1) that this approximation 
can be generalized to the relativistic MHD regime, but we also 
demonstrate (§4.1.1) how any inherent time dependence can be 
recovered. 

The plan of the paper is as follows. In §2 we present the equa- 
tions of time-dependent relativistic MHD, simplify them using 
the frozen-pulse approximation, and discuss what effect each 
of the various forces acting on the plasma has on the flow ac- 
celeration and collimation. In §3 we describe the r self-similar 
model and in §4 we give the results of the numerical integration 
of the model equations. We discuss general implications of this 
work to GRB sources and other relativistic jet sources in §5, 
where we also summarize our conclusions. 

2. THE RELATIVISTIC MHD FORMULATION 
2. 1 . Governing Equations 

The stress-energy tensor of relativistic MHD consists of 
three parts — matter (subscript M), radiation (subscript R) 
and electromagnetic fields (subscript EM): T KU = T£ v + T£ v + 
^em = 0)L2,3). The matter component is given by 

3m" = (pm + Pm/c 2 ) U K U v +P u r) KV , where p M c 2 = p c 2 + Pl) e M 
is the total comoving matter energy density, P M is the par- 
ticle pressure, U u = (70,7V) is the fluid four- velocity, and 
r\ KV = diag(-l 1 1 1) is the metric tensor (assuming a flat space- 
time and Cartesian space coordinates Xj, j = 1,2,3). Here P q 
is the baryon rest-mass density, /?o<?m = Pm/(T~ 1) is the in- 
ternal energy density, with T denoting the polytropic index 
(= 4/3 or 5/3 in the limit of an ultrarelativistic or a nonrela- 
tivistic temperature, respectively), V is the three-velocity mea- 
sured in the frame of the central object, and 7=1 /(l - V 2 /c 2 ) 1//2 
is the Lorentz factor. The radiation component is given by 
7£" = (p R + P R /c 2 ) U K U v + P R ri KV , where p R c 2 and P R are, re- 
spectively, the radiation energy density and pressure in the fluid 
rest frame. Radiation forces are typically most important in re- 
gions that are sufficiently optically thick that one can take the 
local radiation field to be nearly isotropic and set p R c 2 = 3P R . 
We are most interested in the ultrarelativistic case T = 4/3, in 
which the matter and radiation can be treated (under optically 
thick conditions) as a single fluid. Thus, we henceforth write 
P = Pm + Pr, P = Pm + Pr, and p Q e = p e M + p R c 2 . 

Introducing the specific (per baryon mass) relativistic en- 
thalpy £c 2 , where 

<» 



p c* r-ipQC 2 

and including the contribution of the electric (E) and magnetic 
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(B) fields (measured in the central-object frame), the compo- 
nents of the total stress-energy tensor take the form ( j , k = 

1,2,3) 

E 2 +B 2 



pOO , 



8tt 



T°j = T jQ = [£,pon 2 V + 



ExB 

47T 



(2a) 



(2b) 



, EjE k + BjB k ( E 2 +B 2 \ 
T jk = Zpol 2 VjV k - 1 k — 1 + ( P+ -5— ) »/>* , (2c) 



47T \ 87T 

with r 00 , cT°-ixj, and T' k representing the energy density, en- 
ergy flux, and spatial stress contributions, respectively. 
The electromagnetic field obeys Maxwell's equations 

A<jr 

V-B = Q, V-E=—f ) , 



_ D 1 dE 4tt _ 
VxB=- — + —J, 

c at c 



VX£: 



1 dB 



(3) 



c dt ' 

where J v = (J°,J) is the four-current (with J°/c representing 
the charge density). Under the assumption of ideal MHD, the 
comoving electric field is zero, which implies 



V 

E = x B . 

c 



(4) 



The baryon mass conservation equation is (poU u ) jU = 0, or 



d 

— + V-V )( 7 po) + 7PoV-U = 0. 



(5) 



In the absence of a gravitational field or any other external 
force, the equations of motion are T1" = 0. The momentum 
conservation equation is given by the k = 1,2,3 components, 

7 po (| + V ■ v) &V) = -VP+ />E+ J xB . (6) 

The entropy conservation equation (the first law of thermody- 
namics) is obtained by setting U K T K „ = 0, 

l + v.v). + ,(| + v.v)(i; 

which can be rewritten (using poe = P/[T— 1]) as 



:0. 



+ V- V 



:0. 



(7) 



One can carry out a partial integration of equations (3)-(7) 
under the assumptions of axisymmetry [in cylindrical coordi- 
nates (05 ,<f>,z), d/d<j> = 0] and of a zero azimuthal electric field 
(Ec/y = 0) if the flow is time-independent (e.g., Bekenstein & 
Oron 1978; Lovelace et al. 1986). We now show that the equa- 
tions describing a highly relativistic MHD "pulse" that could 
be identified with a GRB event may, in fact, be cast in a steady- 
state form. We start by noting that, with the above assumptions 
(d/d<f>= O,Zi0 = 0), the poloidal component of Faraday's law 
implies that the poloidal magnetic field is time independent. If 
the field is anchored in an accretion disk, then the poloidal field 
near the disk surface would remain quasi steady at least on the 
timescale of the local radial inflow (~ 03/|y ra |). In the case of 
a GRB outflow associated with the emptying up (by accretion 
onto a black hole) of a disk of finite size, the poloidal field may 
be expected to change significantly only on the timescale of the 
burst duration. At the end of this time interval, the information 
that the field has changed starts to propagate with at most the 
speed of light (the actual speed of propagation is the fast mag- 
netosonic speed, which is generally < c in a material medium). 



As we are concerned with highly relativistic outflows, it is rea- 
sonable to expect that the poloidal field associated with the out- 
flowing fluid elements that produce the burst will exhibit neg- 
ligible explicit time dependence (d/dt w 0) over the duration 
of the burst. (Note, however, that the azimuthal component of 
the magnetic field, which is related to the Poynting flux, will be 
time dependent.) 

The solenoidal condition on the magnetic field, V • B = 0, 
implies that there is a poloidal magnetic flux function A(03 ,z), 
defined by 27rA = JJ B p - dS, which satisfies 



B: 



B p + B 



B„ 



VA x 
05 



(8) 



where the subscripts p and (f> denote the poloidal and azimuthal 
components, respectively. Furthermore, equation (4) together 
with the condition E<p = implies V p || B p , from which it fol- 
lows that there are functions ^ A and O (whose coordinate de- 
pendence we discuss below) such that 



V = 



*4 



4njp Q 



S+O3O0, 



v n 



47r7po B F 



(9) 



Denoting the arclength along a poloidal fieldline by £(05, z), 
we change variables from (03, z,t) to (A,£,s), with s = ct-l. 
For any function $ = <£>(A , I , s), we can define the operator 



9$ 9$ d$ 
We now rewrite the MHD equations using 



v = v 5 -w 



d_ 

ds 



Equation (8) becomes 



B: 



V V A x 



05 



d__ d_ 

di~ C d~s 



-+B, 



whereas equation (4) gives 



n 

E = — VA, 

c 



E = 



C50 



-B„ 



Faraday's law implies 



\7 s n x V S A = c 
and the continuity equation yields 



d{E+B^} 



ds 



V s .(4.7PoV) + (c-y,)^^-4. 7 Po^ 



(10) 



(11) 



(12a) 



(12b) 



(12c) 



= 0, (12d) 



where, using equation (9), V, • {Att^pqV) = (V^a x V s A) • </>/05. 

Turning now to the momentum conservation equation, we 
employ equation (11) and the fact that B p is time independent 
(dBp/ds = 0) to write the current density in the form 

9 (B* x V s £- 



J-- 



c 

47T 



V, x B+ 



ds 



-E) 



Decomposing the vectors using the local Cartesian basis 
(n=-= 

we get 



B n 



B rh xVJ-E = B d 



(b-V s £}b + {h-V s l)h 
n + BAn-VsOb, 



= - E+B, 



-En 
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and hence 



J = 



4lT 



d (E+Ba,) - dB s 

OS OS 



The charge density, in turn, can be written as 



c 4tt ^ ^ 4ir 



BF 

V s E-(n-V s £) — 
Os 



Substituting these expressions into equation (6), we obtain 

d_ 

'ds' 



1P0 (V ■ %)(frV)+jp (c-V p )—(frV) = 



4i:ds I VA 



VA 



-VP + —r [V, • (SIVA)] VA + - ' . (12e) 

4irc z 4ir 

Finally, the entropy conservation equation transforms into 



V-V s (P/p r )+(c-V p )- 



ds 



■ 0. 



(12f) 



For a highly relativistic poloidal motion, when j p 



(1-VZ/c 



-1/2, 



i 7 » 1, one can simplify the equations by not- 



ing the following: 1) Due to Lorentz contraction, the observed 
width of the outflow is 7 times smaller than its comoving width, 
d/ds ~ -yd/dl As c-V p w c/2 7 2 and V • V = V p <9/<^, one 
gets (c-V p )9/& - (1 /27)V • V < V • V. Thus, all terms in 
equations (12) containing (c-V p ) are negligible. 2) The term 
4ir^podV P l ds in equation (12d) is of the order of (4-rT pa / ~{)V • 
V^7 and is thus negligible in comparison with the first term on 
the left-hand side of this equation. The arguments above were 
originally given in the context of a purely hydrodynamic (HD) 
flowbyPiranetal. (1993). 3) Using 1 -V p /c w l/2 7 2 , V^/c < 

(1-Vp/c 2 ) 1 / 2 = l/ 7p , and equations (9) and (12b), we infer 

-(E + B^/Bp = (l-V p /c)®n/V p -Vt/V p < (05O/c)(l/2 7 2 ), 
which remains <C 1 throughout the flow in view of the scal- 
ing 7 oc 03 (see § 4.1.2). Thus, all terms in equation (12e) that 
contain (E+B^) (i.e., the first two terms on the right-hand side) 
are a factor <~ 7 smaller than the last term on the right-hand side 
and can be neglected. The same is true in equation (12c), which 
implies dSl w jcd[(E + fi ) /®B p ] ~ 7(9(0 / 7 2 ) < O. 4) The 
pressure-gradient force in the momentum conservation equa- 
tion can also be neglected. It is much smaller than the Lorentz 
force in the transfield direction, consistent with the fact that the 
field is everywhere close to being force-free in that direction 
(especially so in the region near the origin, where thermal ef- 
fects are most important). Along the field, there is a force that is 
~ 7 2 times larger (namely, the inertial force component associ- 
ated with the V ■ V£ term), as in the purely HD case examined 
in Piran et al. (1993). In general, the pressure force is important 
up to the slow magnetosonic point, where, for highly relativis- 
tic temperatures, 7 ~ (3/2) 1 / 2 , and its contribution is negligible 
in the highly relativistic regime. One can therefore replace the 
term -VP by -V S P or even completely neglect the pressure- 
gradient force in the momentum equation without introducing 
a significant error. 

With the above approximations, all the d/ds terms in the 
continuity, momentum, and entropy equations can be elimi- 
nated, and the conservation equations simplify to a steady-state 
form. Although the label s remains attached to the V operator, it 

1 Any value of T other than 4/3 or 5/3 would imply a nonadiabatic evolution 
momentum equations. 



now serves only to identify a given outflowing shell (or pulse). 
The motion remains effectively time independent and can be 
described as a frozen pulse whose internal profile is specified 
through the variable s. As we noted in §1, the frozen-pulse 
approximation was first applied to relativistic HD outflows in 
GRB sources by Piran et al. (1993; see also Piran 1999). We 
have now shown that this approximation can be extended to rel- 
ativistic MHD flows. In the remainder of this paper we pursue 
this effectively steady-state formulation, but we return in §4. 1 . 1 
to consider time-dependent effects in GRB outflows. 

The full set of effectively steady-steady equations can be par- 
tially integrated to yield several field-line constants: 

a) Equations (9) and (12c) yield the field angular velocity, 
which equals the matter angular velocity at the footpoint of the 
fieldline at the midplane of the disk, 

*a Ba, 

(13a) 



03 4-Kjpo 03 

b) The continuity equation (12d) and equation (9) imply that 
the mass-to-magnetic flux ratio has the form 

4irypoV p 



yS> A = q A (A,s) = 



B„ 



(13b) 



c) The (f> component of the momentum equation (12e) yields 
the total (kinetic + magnetic) specific angular momentum, 

03Ba 

L = L(A,s)=^V 4> -—^ . (13c) 

d) Dotting V into the momentum equation (12e) gives the total 
energy-to-mass flux ratio pc 2 , where 

p = p(A,s) = ^-—4 . (13d) 

e) The entropy equation (12f) gives the adiabat 

Q = Q(A,s) = ^. (13e) 
Po 

Equation (13e) is the usual poly tropic relation between density 
and pressure, but in the current application the polytropic in- 
dex is only allowed to take the values 4/3 (if the temperatures 
are relativistic, in which case matter and radiation are treated 
as a single fluid) and 5/3 (if the gas is "cold," in which case 
radiation forces can be neglected). 1 

Two integrals remain to be performed, involving the 
Bernoulli and transfield equations. There are correspondingly 
two unknown functions, which we choose to be the cylindrical 
radius of the fieldline in units of the "light cylinder" radius, 

x = 05ft/c, (I 4 ) 
and the "Alfvenic" Mach number (see Michel 1969) 

M = (-/V./B^irpoO 1 ' 2 = *a(£/4^o) 1/2 . (15) 

We define the Alfven lever arm by m A = (L/pny/ 2 [and cor- 
respondingly x A = U5 A n/c=(LQ/pc 2 ) 1 / 2 ] and use it to scale the 
cylindrical radius of the fieldline by introducing 

G = 03/03a =x/x a . (16) 

To obtain nondimensional variables, we adopt a reference 
length 05o and a reference magnetic field Bq and define 



K 

05 2 



03 2 



05 2 G 2 



A- 



B Q m 2 



(17) 



and hence require the incorporation of heating/cooling terms into the entropy and 
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The expressions for the physical quantities in terms of the 
defined variables and the explicit expressions for the Bernoulli 
and transfield equations are given in Appendix A. Except for 
the s label, which serves to identify a given shell (or pulse), 
these equations are precisely those of steady-state, relativis- 
tic MHD. Solving these equations requires the specification of 
seven constraints, of which four are associated with boundary 
conditions at the source and three are determined by the reg- 
ularity requirement at the singular points (the modified-slow, 
Alfven, and modified-fast points; see Vlahakis et al. 2000). 

2.2. Forces in the Poloidal Plane 

The momentum equation (6) can be written as the sum of the 
following force densities (for simplicity we use hereafter the 
term force) 

fe+fr + fc+fy+fp+fe+f^O, (18) 

where 

i T = -i 2 P«{v-y s ov 

f c = CD7 2 MVi/® 



temperature force 
centrifugal force 



h=-i 2 P oi(v-v s )v- 

f P = -V s P 

f E = (V s -E)E/4ir 

f B = (V s x B)x B/Att 



inertial 
force 



: pressure force 
: electric force 
: magnetic force 



The "gamma" force f G further decomposes into two terms: 
;=fc +fc rfl , with 



A _Pot 
2c 2 



7 



(V-\7 S V 2 )V, f c 



The poloidal part of the f/ force is 

-7 2 Po£(VV s )U p = 7 2 Po£ 



V, 



7 4 M 
" 2c 2 



9V VA 



{v-wl)v. 



di p 



where ip is the angle between the poloidal magnetic field 
and the disk (cos?/; = V m /V p ), and the derivative d/d£ = 
cos tp <9/<905 = [sin 9 cos(i/> + #)/05] dj 89) is taken keeping A (and 
s) constant. The radius of curvature of a poloidal fieldline is 

n = (dip/dey 1 . 

2.2.1. Poloidal Acceleration 
The projection of equation (18) along the poloidal flow is 

7 2 A>£ 



dv 2 7 4 MV„ 2 dV 2 



+7 Po£— cos ip-p & 



2c 2 di 
B d(wB^) 



/2 dt 
di 



-7 z p v;^+ 



(20) 



03 1 ■ " di 4irm di 
The terms on the right-hand side of equation (20) are recog- 
nized as f G || , fry, f c y, fp||, and f fi || , respectively, where a sub- 
script || denotes the component of a vector along the poloidal 
fieldline. The first term on the left-hand side of equation (20) 
is — f/||, whereas the second term is -fc„|| (note that f £ j| = 0). 

2 The strong poloidal magnetic fieldline plays the role of the wire. In the cold limit one has M< 1 and (with xa < 1) near the base of the flow, implying 



The magnetic force component f B || decomposes into the az- 

imuthal magnetic pressure gradient -d (Bl /8tt) /8i and the 

magnetic tension —B 2 ^ cosip/AirtS. These two parts cancel each 
other when B^A ,05) oc 1 /03 2 ; if B^A ,03) decreases faster than 
03" 2 then the gradient of the azimuthal magnetic pressure ex- 
ceeds the magnetic tension, resulting in a positive f B ii . 

In the nonrelativistic regime (V«c, i<1,(«1, f G || = 
fr|| = 0), the pressure force f P \\ dominates up to the slow mag- 
netosonic point, but the bulk of the acceleration is either mag- 
netocentrifugal — corresponding to the f C || term, which can be 
interpreted in the "bead on a wire" picture (e.g., Blandford & 
Payne 1982) 2 — or a consequence of the magnetic pressure- 
gradient force fg|| (which near the surface of the disk can be 
interpreted in the "uncoiling spring" picture; e.g., Uchida & 
Shibata 1985). 3 

The magnetic force generally becomes important also in 
flows where the centrifugal acceleration initially dominates: in 
this case the inertia of the centrifugally accelerated gas am- 
plifies the B<p component, and eventually (beyond the Alfven 
point) fg|| becomes the main driving force. This force continues 
to accelerate the flow beyond the classical fast-magnetosonic 
point (which separates the elliptic and hyperbolic regimes of 
the MHD partial differential equations) 4 and all the way to the 
modified fast-magnetosonic singular point (see Vlahakis et al. 
2000). The modified (and not the classical) fast-magnetosonic 
surface has the property that it is a singular surface for the 
steady MHD equations when one solves simultaneously the 
Bernoulli and transfield equations (e.g., Bogovalov 1997). Only 
when the magnetic field geometry is given (and one solves only 
the Bernoulli equation but not the transfield one) does the singu- 
lar surface correspond to the classical fast-magnetosonic point. 
The modified fast-magnetosonic surface coincides with the lim- 
iting characteristic, the "event horizon" for the propagation of 
fast-magnetosonic waves, since only beyond this surface all the 
points inside a fast Mach cone have larger fast Mach numbers 
than at the origin of the cone. At smaller distances, for a part 
of a given cone, the converse is true: the opening angle of the 
fast Mach cone becomes progressively larger as one advances 
inside that part of the cone; consequently, a small disturbance 
in the super-fast regime can affect the entire flow. The situation 
is similar to that of light propagating close to a Kerr black hole, 
where the ergosphere (which corresponds to the classical fast- 
magnetosonic surface) marks the boundary between the ellip- 
tic and hyperbolic regimes, and where the singular event hori- 
zon (which corresponds to the modified fast-magnetosonic sur- 
face) is located within the hyperbolic regime: only for a spheri- 
cally symmetric (Schwarzschild) black hole are the ergosphere 
and event horizon equivalent — in direct analogy with spheri- 
cally symmetric flows, in which the classical and modified fast- 
magnetosonic points coincide (see Sauty et al. 2002). 

In the case where the outflow attains a highly relativistic 
speed, the centrifugal acceleration cannot play an important 
role. This is because the nonnegligible that would be re- 



that the azimuthal field satisfies —GSB^/LQ 1 — G (1 



-x 2 A )+M 2 ■ 



1 and hence that the f s ii term is negligible and that a near-corotation (V^/USQ »J 1 —M/G ~ 1) 



holds. The small value of M in turn implies a large density and hence a measurable thermal pressure, resulting in a nonnegligible pressure force at the base. 

3 In this picture, the winding-up of the fieldlines by the disk rotation produces a large azimuthal magnetic field component that is antiparallel to in the northern 
hemisphere (and parallel to in the southern hemisphere), and a corresponding outward-directed magnetic pressure gradient -VsiB^/Sn). 

4 At this point, static fast waves with wavevectors parallel to V p in the central object's frame can exist (i.e., eq. [C2] with w = , k \\ V p is satisfied). The classical 
fast-magnetosonic point is equivalently defined by the condition that the poloidal proper speed equals the comoving proper phase speed of a fast-magnetosonic wave 
whose wavevector is parallel to V p (see eq. [C3]). 
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quired in this case would constrain the maximum value of the 
poloidal speed: V 2 = c 2 (l - 1 h 2 )-Vj < c 2 -Vj. Therefore, in 
equation (20), f c y w (and also f G || w 0). The i P \\ force can be 

neglected since fj-H /f P || = (yV p /c) > 1. The two remaining 
terms are f r || (a force with a relativistic origin) and f B \\. The 
expressions for these terms in equation (20) (or, equivalently, 
eq. [13d] for the total energy-to-mass flux ratio) indicate that 
the bulk Lorentz factor can increase in response to the decline 
in either the enthalpy-to-rest-mass ratio £ (the thermal accelera- 
tion case) or the Poynting-to-mass flux ratio oc -C3Z?0 (the mag- 
netic acceleration case) along the flow. When the temperature 
is relativistic, the initial acceleration is dominated by the tem- 
perature force, but after £ drops to <~ 1 the magnetic force takes 
over: this is the likely situation in GRB outflows (see VK01 and 
§4.1). 

When the outflow speed is only mildly relativistic, the mag- 
netocentrifugal force may be important during the initial accel- 
eration phase, especially if the temperatures are nonrelativistic; 
this is the situation envisioned in the "cold" relativistic-MHD 
solutions of Li, Chiueh, & Begelman (1992). It is, however, 
also conceivable that the magnetic pressure-gradient force dom- 
inates from the start, as might be the case if the azimuthal field 
component at the disk surface is large enough; this is the picture 
outlined in the presentation of the relativistic solutions derived 
by Contopoulos (1994; see also Paper II). 

We can obtain an expression for fan as follows. By eliminat- 
ing £7 between equations (13c) and (13d), we obtain a relation 
between USVd, and B3Z?,/,: 



1 



Mi -4) 



whose divergence along the flow implies 

<9(a5V ) _ Mi -4) d(® B 4>) 



dt 



dt 



(21) 



(22) 



Employing the relations for f B \\, f c \\ (see eq. [20]), and equa- 
tions (A1)-(A3), we obtain 



*c\\ 



W m l ^ (i-4)(i ~M 2 -x 2 ) V B p 



dt 



l-M 2 -x 2 



(23) 



V P (-B ) 

The first term on the right-hand side of equation (23) can give 
rise to either acceleration (when V<f, decreases along the flow- 
line) or deceleration (when increases, as in the corotation 
regime at the base of the outflow). This term, together with 
fo^y, can lead to a situation in which V p increases (resp., de- 
creases) and Vcf, decreases (resp., increases) while the Lorentz 
factor remains roughly constant. The second term on the right- 
hand side of equation (23), which is proportional to fg|| , demon- 
strates that the centrifugal force also has a magnetic compo- 
nent and accounts for the Poynting-to-kinetic energy transfer 
that underlies the magnetocentrifugal acceleration process (see 
also Contopoulos & Lovelace 1994 for a related discussion). 
The form of this term makes it clear why the centrifugal force 
exceeds the magnetic force during the initial stage of the accel- 
eration, when the flow is still nonrelativistic (i« 1,M« 1, 
with Bp > \B 4> \,V 4> >V P ). 

The conclusion from the above analysis is that, even though 
centrifugal and thermal effects could dominate initially, the 
magnetic force eventually takes over and is responsible for the 
bulk of the acceleration to high terminal speeds. Li et al. (1992) 
described the efficient conversion of Poynting-to-kinetic energy 
fluxes in relativistic MHD outflows in terms of a "magnetic 



nozzle" (see also Camenzind 1989). The preceding discussion 
makes it clear that, in essence, this mechanism represents the 
ability of a collimated hydromagnetic outflow to continue to un- 
dergo magnetic acceleration all the way up to the modified fast- 
magnetosonic point (which could lie well beyond the classical 
fast point). When interpreted in these terms, it is evident that 
this effect is not inherently relativistic — this conclusion has, 
in fact, been verified explicitly in the case of the nonrelativistic 
self-similar solutions constructed by Vlahakis et al. (2000). 

2.2.2. Collimation 

The projection of equation ( 1 8) in the direction perpendicular 
to the poloidal flow is given by equation (A8). 

The two largest terms are the magnetic and electric force 
components, which almost cancel each other. 




FIG. 1. — Sketch of two meridional fieldlines (solid) and three meridional 
current lines (dashed). The currents satisfy | / |i<| / |2< / I3. Given that 
J p = V s X = j^V,/ X 4>, with / = JJ Jp-dS = \^SB^, the meridional 
current lines represent the loci of constant total poloidal current (/ = const). 
The magnetic and electric forces are shown for the current-carrying (/y < 0, 
left fieldline) and return-current (/y > 0, right fieldline) cases. 

The magnetic force in the meridional plane has two parts: 
\j p xB 4> = V s (C5fi ) and \j^xB p . The first term (which 
usually dominates) has components in both the flow (along 
b = Bp/Bp) and the transfield (along h = E/E) directions. 



contributes to the acceleration 



The b component f B || 
when J± > 0, where the subscript _L denotes the vector com- 
ponent along n. (If a thermally dominated acceleration regime 
exists near the base of the flow, it is in principle possible to 
have J± < there, corresponding to an enthalpy-to-Poynting 
energy transfer mediated by the magnetic force; see Paper II.) 

The n component J -^$- acts to collimate or decollimate the 
flow depending on whether the outflow is, respectively, in the 
current-carrying (7|| < 0) or the return-current (J\\ > 0) regime 
(see Fig. 1). The second term in the decomposition of the 
meridional magnetic force is related to the curvature radius 1Z, 

B 2 

^Jcf, x Bp = n-^ (75 -n • Vsln I V,A/C3 1), and is directed along 
±n for 7^ «s 0. 
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The electric field always points in the n direction, but 
the electric force = (J°/c)En could lie along either +n 
or — h, depending on the sign of the charge density J°/c. 
By employing the curvature radius TZ, one can write f £ = 

(^-n-VM\x\7 s A |). When TZ > (a collimating 
flow), the effect of the curvature term in the expression for the 
electric force is to oppose collimation, but it is possible for the 
other term in this expression to dominate, leading to fs± > 0. 
For highly relativistic motion (7^ >• 1), J° = -j-V s (V x B) = 

T J \\ + V -rh~ h B <t> ■ (V, x V p ) - ±B p ■ (V, x V4) « /„, and 
the electric force has the sign of 7|| . In this case acts to decol- 
limate the flow in the current-carrying regime and to collimate 
it in the return-current regime (see Fig. 1). 

For comparison, note that, when the motion is nonrelativis- 
tic (x <C 1), f£ is negligible relative to i B . In this case, a flow 
in the current-carrying regime is easily collimated (with f B ± 
balancing the inertial force -f/_L = 7 2 po£V / 2 /TZ, resulting in a 
nonnegligible value of TZ). In the return-current regime, colli- 
mation (resp., decollimation) is produced if -J^Bp + J^B^ > 
(resp., < 0); see Okamoto (2001). 

3. THE r SELF-SIMILAR MODEL 
3.1. Model Construction 

To obtain semianalytic solutions of the highly nonlinear sys- 
tem of equations (A6) and (A7), it is necessary to make addi- 
tional assumptions: in particular, we look for a way to effect a 
separation of variables. 

The most complicated expression is the one for B^ (eq. 
[Al]). In view of the importance of the azimuthal field com- 
ponent, which plays a crucial and varied role as part of the 
magnetic pressure-gradient, magnetic tension, and centrifugal 
acceleration terms in the momentum equation, the only realistic 
possibility of deriving exact analytic solutions is to assume that 
the M = const, G = const, and x = const surfaces coincide, i.e., 
that M = M{x) , G = G(x) ,x = x( X ) (Vlahakis 1998). We aim to 
find appropriate forms for the functions of A such that the ex- 
pressions (A6) and (A7) become ordinary differential equations 
(ODEs). From an inspection of the Bernoulli equation (A6) we 
conclude that, in order to separate the variables \ an d A and get 
a single equation that only has a \ dependence, it is necessary 
to assume that the (V S A) 2 term is a product of a function of A 
times a function of \. As A is a function of B3/G(x) (see eq. 
[16]), there must exist functions Ti\ ,Ti.2 such that 



v s (- 



= Hi 



I ) n 2 (G } . 



There always exist the trivial possibilities G oc r in spherical 
coordinates (r,6,<fi) [A = A(0) when the field is radial], and 
G = G(03) (A = A(03) for a cylindrical field), which are not of 
interest here. After some algebra one can prove that the only 
nontrivial case is to have G = G(0), i.e., x = 0. It thus appears 
that, to obtain an analytic adiabatic solution, it is necessary to 
assume r self-similarity. 

The remaining assumptions for constructing an r self-similar 
solution are that the cylindrical distance (in units of the 
Alfvenic lever arm), the poloidal Alfvenic Mach number, and 
the relativistic specific enthalpy are also functions of 9 only: 



x = x{6), M = M(9), £ = £(0) (with the result for £ following 
from the nonlinearity of the expression for M; see eq. [A5]). 

Following the algorithm described in Vlahakis & Tsinganos 
(1998), we change variables from (r ,0) to (a ,0) (see eq. [17]) 
and obtain the forms of the integrals under the assumption of 
separability in a and 0. The results are given in equation (Bl). 5 

Among the five unknown functions of [G(0), ip(9), M(9), 
x(0), and £ (0)] 6 there are three algebraic equations (eqs. [B2a], 
[B2b], and [B2c]) and two first-order ODEs (eqs. [B2d] and 
[B2e]). After solving for these functions, the physical quanti- 
ties can be recovered using 



B 
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(24a) 
(24b) 

(24c) 

(24d) 

(24e) 



where b = B p /B p = zcos-#+Oi5sin?9 = rcos(0-i?)-0sin(0-??) 
is the unit vector along the poloidal fieldline, h = - VA / | VA |= 
z sin d - 03 cos 1? = — f sin (0 - ■&) - 9 cos {9-d) (already defined in 
§2.1) is the unit vector in the transfield direction in the merid- 
ional plane (toward the axis of rotation), and 1? = 7r/2-r/> is 
the opening half-angle of the outflow (the angle between the 
poloidal fieldline and the axis of rotation). 



z=z 




equator 







-> 
CD 



FIG. 2. — Sketch of r self-similar fieldlines in the meridional plane. For 
any two fieldlines Ay and A2, the ratio of cylindrical distances for points 
corresponding to a given value of 9 is the same for all the cones 8 = const: 
Cl/032 = (ai/"2) 1//2 - 



5 The nonrelativistic limit of our model is not the generalization of the Blandford & Payne (1982) model, examined, e.g., in Vlahakis et al. (2000). The 
nonrelativistic limit can, however, be obtained from the analysis of Vlahakis & Tsinganos (1998): it corresponds to the third line of their Table 3 (setting 
xi = F-2.x 2 =F-5/2,E 2 =Ci =D 2 =0, and ignoring gravity, so it is possible to assume a polytropic equation of state). 

6 Note that these quantities also have an s dependence; see § 4.1.1. 
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The r dependence of all the physical quantities can be in- 
ferred from the expressions (24) on the basis of the known r 
dependence of a (oc r 2 ; see eq. [17]). This is a general charac- 
teristic of r self-similar models. 

The parameters of the model are F and T (=4/3 in this 
study), whereas the "constants" x A (s), n(s), a M (s), q(s), and 
Bo(s)U5l~ F (s) together with two "initial" conditions (corre- 
sponding to the two first-order ODEs) are related to seven 
boundary conditions, as we prove in Appendix B. 

The r self-similar character of the poloidal field-line shape is 
shown in Figure 2. 

The model described above is the generalization to a "hot" 
(£ > 1) gas of the only known exact semianalytic solution of 
the relativistic MHD equations, the "cold" r self-similar wind 
solution found independently by Li et al. (1992) and Contopou- 
los (1994). (The force-free model of Contopoulos 1995 can be 
regarded as a special case of the latter solutions, corresponding 
toM = 0andx A = l.) 7 

3.2. Singular Points 

3.2.1. Alfven Singular Point 

It is obvious from equations (24) that the Alfven point, where 
G 2 =M 2 +x 2 = 1, is singular. At this point 



with B being the magnetic field amplitude and with the square 
of the proper sound speed given by 



U 



2 (r-i)(g-i) _ 
- (2-rx+r-i \-c 2 /c 



.2 ' 



(27) 



2 B 2 (l-x 2 ) 
hV e ) 2 = 



In fact, as V„ II B„ 



Singular points appear wherever D = 0; these are the mod- 
ified slow and fast-magnetosonic singular points. They cor- 
respond to points where static slow/fast-magnetosonic waves 
with wavevectors along 9 in the central object's frame can ex- 
ist (i.e., eq. [C2] with u) = 0,fc || 9 is satisfied). 9 The modified 
singular surfaces, which correspond to the "limiting character- 
istics" of the self-similar flow, have previously been considered 
in connection with the nonrelativistic solutions (e.g., Blandford 
& Payne 1982; Tsinganos et al. 1996; Bogovalov 1997; Vla- 
hakis et al. 2000). In order for the solution to pass through a 
singular point where D = 0, N = must simultaneously hold 
(yielding the respective regularity condition). 

3.3. Boundary Conditions and Numerical Integration 

When solving the steady, axisymmetric, ideal MHD equa- 
tions under the assumption that the azimuthal electric field (as 
measured in the central object's frame) vanishes (E^ = 0), seven 

(25) integrations are required, corresponding to seven unknowns: 
the gas density and pressure, the three components of the ve- 
locity, and two functions related to the magnetic field (e.g., A 
and B$). 

Correspondingly, seven boundary conditions determine a 
unique solution. Five of them are the integrals ^ A , fi, L, \i, and 
Q, which, as discussed in §2.1, are conserved quantities along 
the meridional fieldlines. 10 The other two correspond to "initial 
conditions" on the functions G and M, which are obtained by 
integrating equations (A6) and (A7). 

In a physically viable solution, the flow starts with a 
sub-slow-magnetosonic velocity and must satisfy the causal- 
ity principle: any disturbance in the asymptotic regime can- 
not influence the solution near the origin through magne- 
tosonic or Alfven waves. Since the flow starts with a small 
velocity, it must cross three singular surfaces: the modi- 
fied slow-magnetosonic, the Alfven, and the modified fast- 
magnetosonic. 11 The related three regularity conditions are ef- 
fectively three boundary conditions that must be satisfied in 
order for the solution to pass smoothly through the singular 
points. Implementing this procedure is a highly nonlinear task, 
since the positions of the singular points are not known a pri- 
ori and must be obtained simultaneously with the solution. All 
in all, only four boundary conditions remain free and can be 
specified (e.g., on a surface near the origin of the flow). 
Us Bg (1 —x ) i n me r self-similar model that we examine, in which we 

end up integrating ODEs in the variable 9, it is convenient 

(26) to choose as the initial surface a cone 9 = 9 t (where here and 

7 The r self-similarity was first employed by Bardeen & Berger (1978), who examined purely HD flows, but it has become well-known only after Blandford & Payne 
(1982) used it to construct a nonrelativistic MHD disk-wind model. The latter work has subsequently been generalized by many authors (see Vlahakis et al. 2000 and 
references therein). 

8 An equivalent statement is that the Alfven surface marks the locus of points where the flow proper velocity in any direction in the meridional plane is equal to the 
comoving proper phase speed of an Alfven wave that propagates in that direction (see eq. [C3]). 

9 An equivalent statement is that, at these singular points, the magnitude of the flow proper velocity along 9 is equal to the comoving proper phase speed of a 
slow/fast-magnetosonic wave propagating along 8 (see eq. [C3]). 

10 These quantities can be regarded as Riemman invariants; the corresponding characteristics all coincide with the meridional fieldline. 

11 As noted in §2.2.1, the latter surface represents the "event horizon" for the propagation of the fastest waves. The Alfven surface plays a similar role for the Alfven 
waves, but the slow-magnetosonic singular surface is not the "event horizon" for the propagation of slow-magnetosonic waves; it is just the limiting characteristic in 
the sub-slow hyperbolic regime (see Vlahakis 1998). 



47T/90C 

this relation holds not only for the 9- 
components of (V ,B), but for components in any direction in 
the meridional plane. Thus, on the Alfven surface, static Alfven 
waves with wavevector in any direction in the meridional plane 
(in the central object's frame) can exist (i.e., eq. [CI] with u> = 
and = is satisfied). 8 

In order for the solution to pass through the Alfven singular 
point, the Alfven regularity condition (B6) must be satisfied. 
The latter determines the slope of the "Alfvenic" Mach number 
at the Alfven point, (dM 2 /dO) = p A , which is related to the 
Alfvenic value of the magnetization function a A (see Appendix 
B). 

It is seen from equation (25) that the Alfven point is always 
located inside the light surface x = 1 . Note, however, that if 
x A w 1~ (corresponding to the force-free limit, M w 0), the 
Alfven and light surfaces almost coincide. 

3.2.2. Magnetosonic Singular Points 

It is straightforward to obtain an expression for dtp jd9 as a 
function of dM 2 jd9 using the derivative of the Bernoulli equa- 
tion (B2c). After substituting in the transfield equation (B2e), 
the latter takes the form dM 2 jd9 = M jD, where the denomina- 
tor can be written as 



7- 



U} B 2 - 



c 2 4irpot;c 2 ) c 2 4irpot;c 2 



RELATIVISTIC MHD AND GRB OUTFLOWS: I. 



9 



in what follows, a subscript ; denotes an initial value). One 
can start the integration by specifying seven initial conditions 
(i.e., seven functions of r) on this cone (where r is the ar- 
clength along the conical surface; it should not to be confused 
with the distance from the central object, which is given by 
[(rsin#,) 2 + (rcos6','-z c ) 2 ] 1/ ' 2 — see Fig. 2). For example, one 
can specify F , T and Cj > , j = 1 , . . . , 7 such that Bg(r , #,) = 
-d r F ~\ B^r , 6d = -C 2 r F ~ 2 , V r (r , Bd/c = C 3 , V e (r , Od/c = -C 4 , 
V+irM/c = Cs, po(r,6i) = C 6 r 2(F ~ 2 \ and P(r,0 ( ) = C 7 r 2 ^ 
(with the C2, ■ ■ ■ fii being in general functions of s). Note that, 
in the framework of this self-similar model, the specified func- 
tions of r must be consistent with equations (24); if they are 
not, the given scalings will not be reproduced on subsequent 
(6 > 9f) cones. 12 By inverting the system of equations (24), one 
can obtain x A , a M , q, G(6i), M(6>,), ip(0i), and fio03 2 1 ~ F (see eqs. 
[B7]). 13 Three of these parameters are adjusted to satisfy the 
regularity conditions at the three singular points. We recall in 
this connection that T and F are regarded in our formulation as 
model parameters (see §3.1), and we note that ip(6i), Xa, <Jm, 
and q correspond to the fieldline constants /j,, xa, ^a, and Q, 
respectively. 

Next we describe our numerical approach. We have found it 
most convenient to start the integration from the Alfven point, 
since this makes it easier to satisfy the Alfven regularity con- 
dition. We choose a small angular interval < d9 <C 1 and 
specify the model parameters F , T together with the following 
six parameters: 9 A , ipA, xa, <Jm, q and Botsl~ F . (The latter pa- 
rameter does not appear in the system of equations [B2], but 
it affects the magnitudes of the electromagnetic field, density, 
and pressure through eqs. [24].) The seventh parameter is pa, 
which is given from the Alfven regularity condition (see Ap- 
pendix B). We start the integration from 9 = 9 A ± d9, setting 
M 2 = 1 -x\ ± p A d6, G 2 = 1 ± 2cos V'Asin" 1 9 A cos" 1 (tp A + A )d9, 
and ip = ipA- Using the upper (lower) signs, we integrate up- 
stream (downstream) from the Alfven point. Before the first 
step, we evaluate the parameter p, (which is used along the in- 
tegration path to yield iff) from equation (B2c). Whenever the 
solution hits a singular point, we adjust one of the above six 
parameters until a smooth crossing is achieved. 

In our model we have ignored gravity, and consequently we 
expect the flow to be nonsteady in the sub-slow-magnetosonic 
regime. We therefore do not attempt to obtain steady trans- 
slow-magnetosonic solutions, and thus we do not continue the 
integration upstream of the slow-magnetosonic point. This does 
not, however, affect our ability to study magnetic effects, as 
these only become important downstream of this singular point. 

3.3.1. The Roles of Y, F, and z c 

The polytropic index T controls the thermodynamics of the 
flow. For adiabatic flow problems such as the one considered in 
this paper, T = 4/3or5/3 depending on whether the tempera- 
ture is relativistic or not. 

The exponent F controls the current distribution. The 
poloidal current lines are / = |G3Z?0 = const , where, by equa- 
tion (24a), 

1 - G 2 



I = -cW B c 



(28) 



2Fa M l-M 2 -x 2 
Thus, for F > 1, the current | / | is an increasing function of 

12 This is a good potential test for numerical codes solving steady-state equations: 
reproduce the self-similar solution. 

13 One can infer /1 from eq. [B2c] and use it in place of if)(6i) as a fieldline constant. 



a near the base of a trans- Alfvenic flow (where jjjj^s ~ 1)> 
corresponding to the current-carrying regime (see the left field- 
line in Fig. 1). The larger the value of F— 1, the stronger 
the magnetic pinching force that collimates the flow, and hence 
the faster the collimation. The case F < 1 corresponds to the 
return-current regime (represented by the right fieldline in Fig. 
1), whereas F = 1 corresponds to radial meridional current lines. 

In this paper (as well as in Paper II) we concentrate on the 
case F > 1 , which should provide a good representation of the 
conditions near the axis of highly collimated flows such as GRB 
jets. However, in view of the inherent simplifications of the 
self-similar formulation, this choice is not unique. For com- 
parison, we present an F < 1 solution in §4.2.3, and we also 
employ a solution of this type in a forthcoming publication (N. 
Vlahakis & A. Konigl, in preparation) in which we model rel- 
ativistic AGN outflows. A realistic global solution would en- 
compass both the current-carrying and return-current regimes, 
as sketched in Figure 1 . This situation might be mimicked with 
the help of two, properly joined, self-similar models: one (with 
F > 1) that applies near the axis, and the other (with F < 1) that 
applies further out (at larger cylindrical distances). 

The parameter z c (the z coordinate of the disk in the given 
system of coordinates; see Fig. 2) affects only the boundary 
conditions on the disk. For example, z c can be used to mimic a 
Keplerian rotation law. (Recall from § 1 that, in the relativistic 
case considered here, one cannot naturally incorporate a Kep- 
lerian rotation law as in the nonrelativistic r self-similar solu- 
tions.) In our model, il = cx{9)/W, thus, for z c = 0, f2 cx 1/C3 
along the conical surface of the disk. However, for z c > 0, 
points on the surface of the disk at different cylindrical dis- 
tances 03 correspond to different values of 6>(05) (see Fig. 2), 
and ft decreases faster than cx 1/03 (with the rate depending on 
how fast the function x decreases along the surface of the disk). 



4. RESULTS 

The solutions we present in this section are motivated by the 
GRB outflow scenario outlined in § 1 . We approximate the out- 
flow from a disk around a stellar-mass black hole as a pair of 
"pancakes" (see, e.g., Piran 1999) that move in opposite direc- 
tions away from the disk surfaces. The flow originates from the 
inner part of the disk, which extends from the last stable orbit 
around the black hole at r; n to an outer radius r out , which for def- 
initeness we choose to be 2>r m . For simplicity we set z c = 0, so 
(see eq. [17]) a out /a m = 9. We take into account the baryonic 
matter, the electron-positron/photon fluid, and the large-scale 
electromagnetic field. 

As noted in §3.3.1, we focus on solutions in the current- 
carrying regime (F > 1). For this choice of F, equation (28) 
implies that the current vanishes smoothly as the axis (a = 0) is 
approached. A general property of our solutions is that the flow 
reaches an asymptotic cylindrical region where the acceleration 
terminates (see § 4.1). Since we seek to maximize the accel- 
eration efficiency, we consider flows that do not collimate too 
fast (see § 2.2.2), and hence we focus on the smallest possible 
values of F. We therefore choose F > 1. 

Near the origin, the thermal energy associated with the ra- 
diation and e ± pairs is nonnegligible; the optical depth is then 
large enough to ensure local thermodynamic equilibrium. We 

starting with the specified forms of the physical quantities on a cone, they must 



10 



VLAHAKIS AND KONIGL 



therefore assume that the gas (consisting of baryons with their 
neutralizing electrons as well as photons and pairs) evolves adi- 
abatically. The polytropic index is fixed at T = 4/3, correspond- 
ing to relativistic temperatures for the matter and to a blackbody 
distribution for the radiation. Using the Stephan-Boltzmann 
constant 3P R /T 4 = 3P/T 4 (1 +P m /Pr) (where T is the tempera- 
ture) and equations (24e) and (B2b), the temperature in units of 
the electron rest energy is 



Novikov, & Paczyhski 1991) 



e = 



(F-2)/4 



? 1 /4(i + p M /p R )i/4 ( /r C7M )i/2 ^1.25 x 10 13 G 



Bn 



1/2 



(29) 



or, equivalently,^= l+(l+PM/^R)e 4 (po/1-39 x 10 4 g cm" 3 ) 
The matter-to-radiation pressure ratio Pm/Pr is constant in the 
large-temperature limit 8> 1, where, given that the pair num- 
ber density greatly exceeds the baryon number density, the pair 
distribution may be approximated by a Maxwellian with zero 
chemical potential. In this case Pm/Pr = 180/7T 4 w 1.85 (a 
value very close to the more accurate Pm/Pr = 7/4 = 1.75 that 
characterizes a Fermi distribution). In the large-temperature 
limit the pair number density has the T = 4/3 polytropic scal- 
ing (oc 8 3 ), but at < 1 it decreases exponentially, resulting 
in Pm/Pr <C 1. Our polytropic model captures both limits but 
not the intermediate temperatures. Which of the two approxi- 
mations (Pm/Pr = 1-85 or 0) is more accurate depends on the 
initial temperature 9,. For 0, > 1 we choose Pm/Pr = 1-85, 
which yields the correct pressure-temperature relation during 
the initial stage of the flow, when thermal effects are important. 
This choice introduces an error (leading to an underestimate 
of the Lorenz factor) in regions where < 1 but the radia- 
tion energy is nonnegligible in comparison with the baryon rest 
energy. However, because of the weak (a power of 1/4) depen- 
dence of on I + Pm/Pr in equation (29), this error remains 
small. We note in this connection that we also neglect the pair 
rest-energy density, since it is much smaller than the matter 
pressure in the 8> 1 regime where the pair contribution is 
maximized. 

The mass-loss rate in the outflow is M = 2 JJ jpoV ■ dS = 
ff~*AdA, or, 



a; 



F-l 



■a: 



F-l 



M = 



2Fumc 



In 



F-l 

Q-out 



F4\ 



F = l. 



(30) 



dr = 7 



1 COSOJ£ 



n e aj- 



d{± 



(31) 



If Af is the burst duration, then M b ~MAt is the total baryon 
mass ejected. The total energy is ss fjM b c 2 , and initially it 
resides predominantly in the electromagnetic field; the initial 
thermal energy is approximately £;M/,c 2 w (£,-/ /i)£;. 14 

In VK01 we considered two lower limits on the baryon load- 
ing, corresponding, respectively, to the requirements (1) that 
the flow be optically thick in the region where the pair number 
density drops below that of the baryons and (2) that the flow be 
matter-dominated when it becomes optically thin. We also ob- 
tained an upper limit on po from the requirement that the flow 
be optically thin in the internal GRB-shock regime. In the so- 
lutions presented in this paper we evaluate the optical depths 
more accurately. Specifically, consider two neighboring field- 
lines labeled by A and A+dA. Along a direction ( that makes an 
angle luq to the flow velocity, the optical depth is (Abramowicz, 

14 Different shells have different baryon mass densities, so the more accurate expressions are M b 
to-kinetic-energy ratio) / t;j(s)"/j(s)M(s)c 2 ds/c. 



sino^ 

where n e is the electron/positron number density, oj is the 
Thomson cross section, and d(± is the distance between 
the fieldlines. The optical depth is minimized when ui^ = 
arcsin(l/7), for which dr = n e a T d(± (corresponding to pho- 
tons moving perpendicular to the flow direction in the co- 
moving frame). Starting from a point on the inner fieldline, 
we chart the photon trajectory by using u)£ = arcsin(l/7) and 
dr = n e a T d(± until the outer fieldline is reached, and use this 
information to evaluate the optical depth. 

We now present the results of the numerical integration for 
four representative solutions (labeled a, b, c, and d), for which 
the parameters are given in Table 1 . The most important phys- 
ical quantities are plotted in Figure 3, in which each column 
corresponds to a given solution. The properties of these solu- 
tions are described in detail in the following subsections. 

4. 1 . Solution a: A Hot, Fast-Rotator Outflow 

This solution represents our fiducial model of a trans- 
Alfvenic GRB outflow in which the Poynting flux exceeds the 
enthalpy flux at the origin (p 3> £,•). The flow starts from 
the disk with a nonrelativistic velocity and in a short distance 
crosses the slow magnetosonic point (where V p w c/\/3). The 
acceleration in this regime is due to the pressure gradient force 
(the centrifugal acceleration is negligible for £, 3> 1). The slow 
magnetosonic point arises from the interplay between the ver- 
tical gravitational and pressure-gradient forces. As we neglect 
gravity, we start the integration slightly above that slow point 
(7>(3/2) 1/2 ,seeFig. 3a2). 

In this initial sub-Alfvenic regime, Mf < 1 -4,G 2 < 
1,4<4 Equation (A2) for the Lorentz factor gives 

67(l-^)~M(l-4), or (32) 
Thus, a Poynting-dominated flow (p ^> £,) is always close to 
being force-free (x A w 1), and the initial enthalpy-to-Poynting 
flux ratio satisfies w 1 -x\ (= 0.01 in the displayed solu- 
tion). 

As the flow moves downstream it crosses the Alfven singular 
point. We solve numerically for the slope p A of the Alfvenic 
Mach number that satisfies the Alfven regularity condition (see 
§3.2.1). 

In the super-Alfvenic regime there are three possible cases: 

1. The flow recollimates (tp > ir/2) and there is a termination 
point where the entire energy is suddenly transformed into ki- 
netic motion (corresponding to the solution hitting the modified 
fast singular surface but not being able to cross it); 

2. The flow starts to decelerate at some heigh above the disk; 

3. The flow becomes asymptotically cylindrical. 
The last case is the only physically acceptable one and corre- 
sponds to a specific value of one of the model parameters. (For 
example, for the chosen parameter set F, A , p A , x A , <jm, and q, 
there is a unique value of p that corresponds to the cylindrical 
solution). This is the "magnetic nozzle" mechanism first de- 
scribed by Li et al. (1992) and discussed in §2.2.1. When the 
flow has cylindrical asymptotics, the asymptotic regime 9 = 
is the only possible solution of T> = (where T> is given by 
eq. [26]). In this case, the modified fast-magnetosonic singular 

J M(s) ds/c, £j = J [i(s)M(s)c 2 ds/c, and (for the initial enthalpy- 
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solution a solution b solution c solution d 

G3(cm) d(cm) 03(cm) CO(cm) 

10 s 10 7 10 B 10' 10'° 10"l0 6 10 7 10 8 10 s 10'° io"io" 10' 10" 10 s 10'° io'io" 10 7 10" 10 9 10'° 10 




FIG. 3. — Main properties of the four representative solutions discussed in §4. See text for details. 
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Table 1 

Parameters of Representative Solutions 1 " 



solution 


F 


4 


cm 


1 


M 


floral (cgs) 


0a(°) 




PA 


&A 


a3o(10" cm) 




1.01 


0.99 


5000 


37.396 


10116.1 


1.35 x 10 2 " 


35 


72.7 


-3.9 x 1Q- 1 


86.6 


2.5 


b 


1.01 


0.9999 


5000 





9997.4 


1.05 x 10 20 


35 


72.7 


-5.0 x 10~ 4 


6820 


2.5 


c 


1.01 


0.15 


300 


322501 


4053.2 


5.93 x 10 19 


35 


72.0 


-3.0 


0.18 


3 


d 


0.7 


0.96 


1000 


44.92 


2150 


6.83 x 10 21 


10 


83.4 


-0.93 


23 


2.5 



T In all cases f = 4/3 , Zc = 0, and /-,,„,/>;„ = 3. 
* The solution presented in VK01 is the same as solution a, except that B ®l~ F = 2.96 X 10 19 cgs and Ho = 7.8125 X 10 5 cm. 



surface (the "event horizon" for the propagation of fast waves) 
corresponds to the asymptotic cylindrical regime (where all the 
fast Mach cones along each fieldline have the same opening an- 
gle). As in the "cold" solutions derived by Li et al. (1992), the 
critical value of p, is always close to 2a M - 

Figure 3al shows the meridional projections of the innermost 
and outermost fieldlines on a logarithmic scale. Our chosen 
initial cylindrical distances are 03, « 1.7 x 10 6 cm for the in- 
nermost fieldline and s=a 5.2 x 10 6 cm for the outermost one; 
these distances are only slightly larger than those of the foot- 
points of these fieldlines in the disk. Figure 3al also depicts 
the optical paths of photons that originate at three points on the 
innermost fieldline; these were calculated according to the pro- 
cedure described at the beginning of this section and are plot- 
ted as dashed lines labeled 1, 2, and 3. Figure 3a7 shows the 
corresponding optical depths t(1), t(2), and t(3) along these 
paths as a function of r out /r m = (aw/c^in) 1 • It is seen that, as 
one moves downstream (with the initial point on the innermost 
fieldline shifting to a larger height z above the disk), the optical 
depth gets smaller, and the flow eventually becomes optically 
thin (r(3) w 1 for r out /r in = 3). 



10" 



10" 



10~' 



! along V p j 


along -VA 

\ ' \s I 


r <c\\ X 

n \ X 




- _f /wb+'e ' 

r v\ ■ 1 

v\ ''■ 

\V\ 1 


\ \\ 

r \n \\ i 

[ \ \\ 

! Y* \\ 

- \M v 3 


\ \ 
\\\ 

\\\ 

- -fp\\ 1 

\ \ 


; % ; 


-X 




: (b> j 



10' 10" 



rate 



10'" 10' 



10" 



C5(cniP 



10" 



FIG. 4. — Force densities in the meridional plane (a) along the poloidal flow 
and (b) in the transfield direction (toward the axis) for solution a. 

Figure 3a7 also shows the radial profile of the mass-loss rate. 
For r out /r m = 3, M rs 10" 7 M Q s" 1 , corresponding to a total 
ejected baryonic mass of M/, m 10" 6 M Q for a typical burst du- 
ration Af m 10 s. 

The total energy flux along the poloidal flow is pc 2 x -fpoV p 
(where -fpoV p is the mass flux). The main part of the total en- 
ergy flux (at least in the initial phase of the flow) is the Poynting 



flux (see eq. [2b]) 
V 



(Note that this term is positive since < 0.) The remaining 
part (p^ + WV.B^/^a) x -fpoVp = £-fc 2 x jpoVp (see the energy 
conservation relation, eq. [13d]) includes the kinetic energy 
flux (7Poc 2 (7~ l)Vp) an d the enthalpy flux. The total energy 
loss rate is £i = 2ff jp V p (£jc 2 - OJ^B^/^a) ^-dS = pMc 2 , 
and for Af w 10s the total energy injected into the two oppo- 
sitely directed jets is w 1.8 x 10 52 ergs. 

Figure 3a2 shows the various energy fluxes in units of 
7Poc 2 Vp (the mass flux xc 2 ) as well as the Lorentz factor 7 
as functions of 03, the distance from the axis of rotation along 
the innermost fieldline. We distinguish three different regimes: 

1 . Thermal acceleration region: 

From the initial point (slightly above the slow point) up to ~ 10 8 
cm, £7 « £ ; = const . (The exact initial values for solution a are 
7; = 1 .2 and = 1 14.8.) The acceleration is primarily thermal: 
enthalpy is transformed into kinetic energy. In Figure 4a we see 
that the force — fgy (which measures the increase in 7) is equal 
to the temperature force f r || (which describes the decrease of 
£). The Poynting-to-mass flux ratio (-ESflB $ / a) is essentially 
constant, which means that the field is force-free; it only guides 
the flow. To a good approximation, 7 sd 7,03/03, and £ w £,03, /G3. 
So long as £ 3> 1, equation (29) gives £ oc 9; thus, 9 w 9,03,703. 
This is verified in Figure 3a3, which shows the poloidal mag- 
netic field (in units of 10 14 G) and the baryon rest-mass den- 
sity (in units of 100 g cm" 3 ). It is seen that po oc 1/03 3 (as 

expected from the poly tropic relation P oc p 4 / 3 ) and B p oc 1/03 2 
(as expected from the constancy of the mass-to-magnetic flux 
ratio 4-n-fpoVp/Bp = ^ A ). As the azimuthal velocity is negli- 
gible (for a highly relativistic poloidal motion), equation (13a) 
implies O3Z?0 = const., as expected in the force-free regime. In 
general, = cx + VpB^/Bp, so —B$ fis xB p = E, consistent with 
Bp oc 1 /03 2 and B^ oc 1 /03. 

Figure 3a4 shows the angle between the total magnetic field 
and its azimuthal component (which, given that B p oc 1/05 2 and 
-B^ oc 1/03, decreases as 1 /03), the "causal connection" open- 
ing angle arcsin(l/7) oc 1/03, and the opening half-angle of the 
outflow i3 = 7r/2 — if} (its initial value is i?,- fia 25°). 
Figure 3a5 shows the pressures associated with the poloidal, 
azimuthal, and comoving magnetic field components, as well 
as the total thermal (lepton + radiation) pressure contribution 
(B 2 /8tt, Bl/in — dashed lines; B 2 /8tt = (B 2 -£ 2 )/8tt,P — 
solid lines, respectively). In general, B 2 Q fis B 2 +B 2 b /'j 2 fia 
B 2 (l+x 2 /7 2 ). As x and 7 are both proportional to 03, their 



XSVLBa 



ratio is a constant, x/~/ = JCj/7; 
verified by the figure. 



' Xj <C 1. Thus, B co w B„, as 
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Figure 3a6 shows x 2 = (mfl/c) and M 2 = (jV p ) / (B 2 p /4np () 
(the squares of the "light cylinder" radius and the Alfvenic 
Mach number, respectively) as well as the square of the "classi- 
cal fast-magnetosonic proper Mach number" M 2 = (jV p 
where Uf is the larger solution of the quadratic 



> 2 



/' 



U, 



Ui 



u 



B 1 



U? B 2 ,(1- X 2 ) 
c 2 ' 4irpo£c 2 J ' c 2 4-Trpo^c 2 



= 0. 

(34) 

The point where x = 1 corresponds to the light surface, which is 
close to the Alfven surface x = xa (more accurately, the Alfven 
surface is where M 2 = 1 — xfy. 

The point Mf = 1 is the classical fast-magnetosonic point. It is 
seen that, up to that point, M w Mf. This can be understood 
by noting that, since B 2 -E 2 w B 2 > 4tt pq£U 2 , the solution of 
equation (34) yields M 2 /M 2 f = l/2 + (l/4+x 2 M 2 U 2 /j 2 V 2 ) 1 / 2 . 
It follows that, so long as x 2 M 2 U 2 /^ 2 V 2 < x 2 M 2 /2 < 1 /2 (us- 
ing U 2 < c 2 /2), Mf w M. 

At the classical fast point (subscript /) M w 1 <C x, so equation 
(A2) for the Lorentz factor gives 7/ w p/^jx 2 . If this point is 
located inside the thermally dominated region (£/ 3> 1, as in 
the depicted solution), then, using 7/£/ w £,7,-, one infers w 
(M/67i) 1/2 - As 7 /x w 7,/xi, it follows that 7/ w (pji/^x 2 )^ 2 . 

2. Magnetic acceleration region: 
From the end of the thermal acceleration zone, where £ w 1, 
up to 03 w 10 1() cm, it is seen from Figure 3a2 that the Lorentz 
factor continues to increase as 7 oc 05^, with (3 a constant w 1. 
This exponent is in general different for different solutions. We 
find that for solutions with a given p but different (3 is al- 
ways less than 1 and decreases with increasing For larger 
£i the magnetic effects are less important, resulting in a weaker 
collimation of the flow and hence in a larger asymptotic width 
03 oo. Since the final Lorentz factor is <~ p/2, the exponent 
(3 = d\n~f / d\nTB should be smaller. Centrifugal acceleration can 
also influence the magnitude of this exponent: for mildly rela- 
tivistic flows we expect this acceleration to increase the lever 
arm of the flow, resulting in a lower value of (3. 
The acceleration in this region is due to magnetic effects: 
Poynting energy is transformed into kinetic energy. Figure 4a 
shows that the force -f G || (which describes the increase in 7) is 
equal to the Lorentz force f B ii (which derives from the decrease 
of I GBcf, |). The Poynting-to-mass flux ratio declines from its 
initial value w pc 2 as pc 2 - | TBSIB^/^a |oc 03^, with (3 w 1, a 
result verified by Figure 3a2. However, the deviation from the 
initial value of pc 2 is not too strong, so the scalings B p oc 1/05 2 , 
po oc 1/03 3 , and -B p /B,p oc 1/05 remain approximately the same 
as at smaller radii (see Figs. 3a3, 3a4, and 3a5). 
Figure 3a6 shows that x S> M but that M increases 
faster than x. As the Poynting-to-matter energy flux ra- 
tio is \c(E x B) ■ b/4n\ / [^ 2 p Q c 2 V p ] = (p - £ 7 )/£ 7 = (x 2 - 



tion (24c), to obtain the fieldline slope, 
Fa M x 2 + M 2 . sin(0-t?) $ 



1- 



1- 



A=const 



p x*- sin 9 

(35) 

(The same result can be found using eq. [B2c].) So long as 
x ^> M, the slope is \—Fgm/p ~ 1/2 (as the critical value for 
p is close to 2<r M , and F w 1). Thus, the shape of the fieldlines 
is parabolic, z oc 05 2 , as Figure 3al verifies (cf. Contopoulos 
1995). 

3. Asymptotic cylindrical region 
At the end of the magnetic acceleration region the flow becomes 
cylindrical: Figure 3al shows that the fieldlines tend to a con- 
stant 05, whereas Figure 3a4 indicates that the opening half- 
angle = 7r/2 — ip tends to zero. 

Figure 4b shows that, although the electric force almost cancels 
the magnetic force (fs± + (e± <C fs±), their sum fg± + (e± ~ 
-fcj_. On the other hand, -f/j_ <C -fc_L, or (using eqs. [A8d] 
and [A8e]), 05/7?. <C V% smip/V 2 , which means that the poloidal 
curvature radius vanishes (i.e., the poloidal fieldlines become 
straight) — a characteristic of cylindrical collimation. (In the 
initial acceleration region near the disk it is seen from Fig. 4b 
that { b ± + (e± ~ -f/_i_, which implies that the fieldlines curve 
toward the axis.) 15 

In the cylindrical region, the transfield force-balance equation 
becomes 

V 2 d / p+ tfV_B|_. 

87T / 47T05 



^0^7 



d_ 

05 dXS 

or, using a oc 05 2 and equations (24), 
,, _2-F 



(36) 




F-\ 
47T05 



(37) 



As the left-hand side of equation (37) is small but positive, 



1) 



(F- 

B 2 p + Bl 



jB\ w + . Using B 2 C0 w B 2 p +B\ / 7 2 and B 2 co = 

E 2 , we conclude that B 2 ^ is invariably > E 2 ; hence, 
only a current-carrying jet (F > 1) can have cylindrical asymp- 
totic s. 

Equation (35) with dtS = gives the final kinetic-to-Poynting 
flux ratio, (M 2 /x 2 ) oc w (p/F<r M )- lwl, corresponding to 
7oc ~ P~Fa M ~ p/2. The value of M 2 becomes as large as 
that of x 2 (see Fig. 3a6), so a rough equipartition holds between 
the kinetic and Poynting fluxes: 700 w (-OSfiB^/^c 2 ) (see 
Fig. 3a2). 

The implied conversion efficiency of <~ 50% between the total 
energy (injected largely in the form of a Poynting flux) and the 
final kinetic energy of the flow is consistent with the internal 
shock scenario for GRBs. Specifically, the asymptotic kinetic 
energy in each jet is ~ 2 x 10 51 ergs for our fiducial parameters. 
With a radiative efficiency of > 10% (e.g., Kobayashi, Piran, 
& Sari 1997; Nakar & Piran 2002b), this implies a radiated 7- 
ray energy of a few times 10 50 ergs, as inferred observationally 
(Frail etal. 2001). 

Figure 3al shows that, as the cylindrical region is approached, 
the Lorentz factor increases more slowly than oc 05, resulting in 
a nonnegligible value of x/7. Since Z? 2 w B 2 (l +x 2 /7 2 ), this 



x 2 i )/(M 2 + 1 —x%) ~x 2 /M 2 , this is another manifestation of the 
Poynting-to-kinetic flux conversion. 

The Bernoulli equation simplifies in this region to 7 w jV p /c, 
which can be used, along with 7 w pM 2 / £(x 2 + M 2 ) and equa- 

15 The acceleration could in principle continue if a transition from positive to negative poloidal curvature were possible. Although the geometry of the r self-similarity 
does not allow such a transition, other types of self-similar solutions can be constructed in which such a transition takes place and the flow continues to accelerate 
(see Vlahakis 2003). 



explains the divergence of the (B 2 -E 2 )/8tt and B 2 /8tt curves 
in Fig. 3a5 near the cylindrical region. (Eq. [37] can be written 
as (F- l)(B 2 -E 2 ) = B 2 +4irp()£j 2 V^+%Tr(2-F)P, which shows 
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that B 2 -E 2 > B 2 p /(F- 1) in the cylindrical regime.) 
The condition that the GRB emission region be optically thin 
to the prompt 7-ray photons typically implies a lower limit 
> 100 on the Lorentz factor of the outflow (e.g., Lithwick & 
Sari 2001). The asymptotic Lorentz factor attained in our fidu- 
cial solution (and, in fact, in the other three solutions presented 
in this section) readily satisfies this requirement. 

4.1.1. Time-Dependent Effects 

As noted at the beginning of §4, we approximate the outflow 
as a pair of pancakes that move in opposite directions away 
from the disk. In the context of the internal-shock scenario, 
each pancake consists of a number N = 100 N2 of shells. The 
width of each shell is 5s = A£/N = cAt/N = 3x 10 9 (Af)i/JV 2 
cm, where the total duration of the burst is At = 10 (Af)i s. 

If the ejection of the shells from the disk surface (£ = 0) starts 
at time t = and ends at time Af , then each shell can be labeled 
by its ejection time f, = s/c, with < s < cAt. The shell "s" 
moves along the poloidal fieldline and at time t its position is 
I = f s / c V p dt. The distance between two neighboring shells s, 
s + Ss is (using V p ~ c-c/27 2 ) 

5£ = 5^ Vpdt^j k-Ss + J^JL^ldt. (38) 

At early times the second term on the right-hand side of equa- 
tion (38) is negligible and the frozen-pulse approximation holds 
(the distance between the two shells is constant, 51 « -5s). 
However, even if the integrand in the second term is tiny, there 
will eventually come a time (call it t c ) when the integral will 
cancel the -5s term (for 5'j/5s> 0, corresponding to the trail- 
ing shell moving faster than the leading one). One can integrate 
the equation of motion for each shell, dz/dt = c(l - I/7 2 ) 1 / 2 , 
with 7 w 7,-BJ/BJ,- ~ Jiiz/zi) 1 ^ 2 , to show that two neighboring 
shells with 5-y <~ 7 will not collide for as long as they move in 
the flow acceleration zone. This is a reflection of the fact that 
the difference in the initial Lorentz factors of the two shells is 
not large enough to compensate for the longer acceleration time 
of the leading shell as it crosses this region, which also implies 
that the frozen-pulse assumption continues to hold throughout 
the acceleration zone. Only after the shells reach the constant- 
velocity, cylindrical-flow region, will the two shells collide 
(5£ = 0). Using equation (38), this will happen at a distance 
i=a "floSs w 2 x 10 16 (Af)i/A^2 cm beyond the end of the acceler- 
ation zone. 

So long as the frozen-pulse approximation is valid, £ = 
f' s / c V p dt w ct-s. As we proved in §2.1, it is possi- 
ble to examine the motion of each shell in this regime 
using the r self-similar model. If we focus on the 
shell so, then the parameter F and boundary conditions 
{x A (s Q ) , fi(s ) , <t m (so) , q(s ) , B (so)Wo(s ) 2 ~ F } that we used to 
specify the solution refer to this particular so, and the cor- 
responding solution {M(8 7 so),x(0 7 so)} is valid only for this 
shell. 

In general, we may choose a different set 
{F ,x A (s), ii(s),cr M (s),q(s),Bo(s)® Q (s) 2 ~ F } for a different 
shell s, but we have to be careful to also satisfy the assumption 
of a quasi-steady poloidal field, which requires all the shells 
to experience the same poloidal magnetic field at any given 
location. This requirement constrains the choice of boundary 
conditions. 



No constraint is necessary in the force-free regime, in which 
the electromagnetic field is effectively decoupled from the mat- 
ter in that there is no feedback from the matter acceleration to 
the field. The electromagnetic field only guides the flow and 
the motion is effectively HD. As we proved in VK01, one sim- 
ply has to replace the spherical r in the radial-outflow scalings 
(e.g., Piran 1999) with the cylindrical radius G5 to get the correct 
scalings in the magnetically guided case. 

However, in the superfast regime (M 2 3> 1 — x\, x 2 ^ x 2 A ), 
one finds that the (appropriately simplified) Bernoulli and trans- 
field force-balance equations (eqs. [B2c] and [B2e]) become 
completely s-independent if one writes M 2 (9,s) = M 2 (9)g(s), 
x 2 (9,s) = X 2 {6)g{s), and n(s)/a M (s) = ^(s )/a M (s Q ). (F 
must also be s-independent.) Thus, we are free to 
specify the functions g(s), x A (s), ctm(^), and q(s) when 
/z(s) = /j,(so)<tm(s) I &m{sq), and from the fact that A is s- 
independent we then also have Bo(s)®l~ F (s)x F (s)/g F / 2 (s) = 
Bo(sq)G5 2 i ~ f (sq)x f ^(sq)/ g F / 2 (sq). These functions correspond to 
the initial conditions for each shell, obeying the quasi-steady 
poloidal magnetic field assumption. 

Using s = ct - £ in any quantity $ = $(A , £ , s), we may find 
either the time dependence following the motion of a particu- 
lar shell: $ = <£>(A , ct - s , s) with s held constant, or the time 
dependence at a given point in space as different shells pass 
by: $ = <&(A,£,ct-£), with £ held fixed. The "initial" value 

= $(A , , s) = $(A , , cti) corresponds to the time t = ti = s/c 
when the shell s is ejected from the surface £ = 0. Thus, we see 
that the s-dependence in $(A ,£,s) represents the initial condi- 
tions at the ejection surface for each shell. 
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FIG. 5. — (a) Baryon mass density plotted as a function of the arclength 
along the innermost fieldline of solution a for several values of time since the 
start of the burst. The pancake width (A£) remains constant cAt. The 
dashed line represents the "time independent" solution for the reference shell 
so = 5 X 10 10 cm, corresponding to the heavy dots, (b) Assumed form of the 
function g(.s)/<r M (.v) across the pancake. This function represents the initial 
baryon mass density in each shell (-fpoi(s) oc g(s) / (Tm(s)\ see eq. [40]). The 
shell s = 0, which is ejected first, has a tiny mass, the following shells have 
progressively larger masses, and the last shells (s < cAt) again have negligi- 
ble masses. 



To illustrate how the time dependence in our model can be 
recovered, we now consider the baryon mass density along the 
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flow. From equation (24d), 

[B (s Q )®t F (s )] * p(s )x 4 A (s )g(s) A 2 ~t 



JPo(r,9,s)-- 



4TTc 2 F7a M (so)g 2 (so) [M 2 (6) + X 2 {0)] a M (s) 



M 2 {0)g{s)-\ +x 2 (s) X 2 (0)g(s) + M 2 (6)g(s) 

M 2 (9)g(s) X\6)g{s) + M\6)g{s)-\- 

Neglecting the second line (which is approximately equal to 1 
in the non-force-free regime), 

g(s)/a M (s) 



JPo(r,0,s) = -fp o (r,8,s Q )- 



(40) 



g(so)/o- M (so) ' 

Figure 5 shows jpo as a function of the arclength along the in- 
ner fieldline at various times (as in the HD case illustrated in 
Fig. la of Piran et al. 1993). At each instant of time t , one can 
identify (using s = ct- I) which shell s is at the point i. Know- 
ing the distribution of 7/50 as a function of I for the reference 
shell sq, one can obtain the density for all the other shells (corre- 
sponding to a particular choice of the function g(s) / ctm(sX such 
as the one depicted in Fig. 5b). 

4. 1 .2. Validity of Assumptions 

The ideal-MHD theory is not always a good approximation; 
in some cases it must be replaced by an exact multi-fluid the- 
ory including a relativistic Ohm's law (e.g., Melatos & Mel- 
rose 1996). A necessary condition for its applicability, derived 
from the need to attain the requisite charge density and cur- 
rent, implies a lower limit on the matter density (e.g., Spruit, 
Daigne, & Drenkhahn 2001). The limits are J <C poq e c/m p , 
J°/c <C poq e /m p (where q e the electron charge), and we have 
verified that they are well satisfied in all of our solutions. In 
particular, the protons and the neutralizing electrons have a suf- 
ficiently large number density to screen the electric field paral- 
lel to the flow and to provide the necessary charge density and 
current. 

Our neglect of gravity has turned out to be a valid approxima- 
tion, since the flow near the disk is thermally dominated and the 
temperature force density near the origin of our fiducial solu- 
tion (at 05, ~ 2 x 10 6 cm, slightly above the slow-magnetosonic 
point) is fj|| Ri 10 19 dyne cm" 3 , which is much larger than the 
gravitational force density associated with a central object of a 
few solar masses. 

In § 4.1.1 we demonstrated that the frozen-pulse assumption 
is valid throughout the acceleration region. 

We have also neglected the terms d(E + B^) /ds and d(E 2 - 
B^/ds in equation (12e). Our results by and large verify 
these approximations, as we find that only near the asymp- 
totic cylindrical region does B co = (B 2 + B 2 ^ - E 2 ) 1 ! 2 deviate 
from Bp (see Fig. 3a5). Even so, the ratio of the terms in 
equation (12e) that contain E+B^ over the last term in that 
equation is j(E + B^/B^, which (using eq. [4]) becomes 
7(1 -Vp/c) < 1. An alternative argument supporting the con- 
clusion that E Ri is based on the requirements that B 2 a > 



andV^ > 0. The former gives (using B 2 = B^ + B^-E 2 and eq. 
[12b]) B^/E 2 > 1 - 1/x 2 , whereas the latter yields (using eq. 
[9]) B\jE 2 < c 2 /V 2 . Thus, the ratio B^/E 2 is bounded in a 
tiny interval around 1, from 1 - 1/x 2 to c 2 /V 2 « 1 + 1 /j 2 . 

4.2. Other Solutions 
4.2.1. Solution b: The Cold Case w 1) 



In the limit of low temperatures (q — > 0, £ — ► 1), we reproduce 
the exact cold relativistic solutions of Li et al. (1992), using a 
different parameter regime appropriate to GRB outflows. Con- 
topoulos (1994) employed the same model but examined more 
massive flows that were not close to being force-free and thus 
were not accelerated as efficiently. 

We present a cold solution in the second column of Figure 3; 
it can be readily compared with the fiducial solution a. 

The cold solution is closer to being force-free, since 7, w 
1 /(l -x 2 ) 1 / 2 (using Vpi w , y 01 - w G3 ( Q), and hence (by eq. [32]) 
p(l -x\) w (1 -x 2 ) 1 / 2 Ri 1. (For comparison, p(l -x 2 A ) rs £ for 
solution a.) 

Since the flow does not pass through a slow magnetosonic 
point in this case, we are able to find the solution all the way 
down to the disk surface. Near the origin the flow corotates 
with the disk (V^ ~ cx) and is accelerated due to the centrifugal 
force f C ||. The Lorentz factor does not increase linearly with 
03 (but rather, to a good approximation, as 7 w 1/(1 -x 2 ) 1 / 2 « 
(1+x 2 ) 1 / 2 ). However, the Lorentz force soon takes over and 
thereafter 7 cx 03 (see Fig. 3b2). 

At the classical fast point Mf Ri M Ri 1, equation (A2) gives 
jf « /i/x 2 , which (using 7/ rj Xf) can be transformed into the 

familiar form 7/ Ri /i 1 / 3 (e.g., Camenzind 1986). Unlike the 
purely monopole case (Michel 1969), this point is located at 
a finite distance from the origin, and the bulk of the magnetic 
acceleration occurs further downstream. 

The flow near the cylindrical regime is similar to that in the 
fiducial solution, as thermal effects are not important anymore. 
However, a signature of the cold solution can be seen in Figure 
3b6 in the weaker deviation of M from Mf. (The solution of eq. 
[34] in this case is M 2 /M 2 = B 2 /(B 2 -E 2 ), and the deviation 
of M from Mf is due only to B 2 -E 2 becoming larger than B 2 , 
which happens near the cylindrical region.) 

4.2.2. Solution c: The Slow-Rotator Case (£,■ Ri p) 

This solution describes a slow rotator, in which the Poynting 
flux is smaller than the enthalpy flux (so p Ri £,). The angu- 
lar velocity of the inner fieldline is il = cxa/GSq = 3873 rad s -1 
(see Table 1), much smaller than the typical value (~ 10 4 rad 
s" 1 ) near a solar-mass black hole. This solution is character- 
ized by significantly higher B p j \ B^ | ratios (see Fig. 3c4) and 
much lower values of x (see Fig. 3c6) in comparison with solu- 
tions a and b. It is also seen that this solution is not force-free: 
x\ = 0.15 <§; 1 (see eq. [32]). The acceleration in this case is 
thermal (£7 Ri p; see Fig. 3c2). 

Nevertheless, since the flow is trans-Alfvenic, the poloidal 
magnetic pressure is much larger than the thermal pressure (see 
Fig. 3c5). This can be understood from the fact that, in general, 
one has at the Alfven point 

(Bl_\ _( 2 (jVp/c) 2 \ ^ 2 7 2 



l-l/£ 1-x 2 j 1-x 2 



> 1. 



(41) 



A \ /A 

Only when the flow is super-Alfvenic everywhere (even near 
the origin) can one obtain a hot, radial, HD-like solution with 

P > B 2 /8tt. 

The magnetic field only guides the flow in this case. The 
collimation, however, is not very efficient (see Fig. 3cl). The 
flowlines "attempt" to become conical, but since in the frame- 
work of r self-similarity the radial distance is given by r = 
G5oa l / 2 (A)G(9)/ sin 6*, the only possible asymptotes as r — > 00 
are — > 6*00. For ^ = we have cylindrical asymptotics (as 
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in solutions a and b), and the only other possibility is for all 
the fiowlines to have the same asymptotic opening half-angle, 
which is the situation depicted in Figure 3c 1 . This behavior can 
be understood if we write the transfield force-balance equation 

as 



b: 



4ttK 



2 



Bi 



+ 



( 



MV, 



47TC5 V V p 



n-ib-n-V s 




(42) 



(see Vlahakis 2003). By neglecting the centrifugal (oc V|) and 
the poloidal magnetic pressure terms and using equation (4) 
(which implies E -B^Vp/c), equation (42) can be rewritten 

as 

4§(* ! « ! - 1 )=-i"^-( 5 r) 2 -^-< 43 > 

The first term on the right-hand side is positive and acts to colli- 
mate the flow (resulting in a positive poloidal curvature, 1Z > 0), 
whereas the second term is negative and has the opposite effect 
on 1Z. In Poynting flux-dominated solutions, the pressure term 
is negligible and the collimation continues (albeit at a very low 
rate; see paper II) until the shape becomes cylindrical. How- 
ever, when the thermal pressure becomes comparable to the co- 
moving magnetic pressure (as in the present solution; see Fig. 
3c5), the two terms on the right-hand side of equation (43) can- 
cel each other. The curvature 1/7?. thus vanishes before the 
cylindrical regime is reached, resulting in an asymptotic coni- 
cal flow. 

4.2.3. Solution d: The Return-Current Regime 

For completeness, we also present a return-current (Jii > 0) 
solution. Since, near the disk surface, -OJfi^ oc a F ~ l is a de- 
creasing function of a in the F < 1 case under consideration, 
the magnetic force acts to t/ecollimate the flow (see § 2.2.2). 
This results in a weaker collimation than in the current-carrying 
(F > 1) solutions (see Fig. 3dl). The electric force succeeds 
in collimating the flow despite the counter effect of the mag- 
netic force; as discussed in §2.2.2, this is associated with the 
positive charge density (J° > 0) expected in the return-current 
regime. Examining the forces in the transfield direction, we find 
that f £ j_ ^> f E ± +{b± ~ _ f/_L ^> -fc_i_; i- e -, the Lorentz force is 
balanced by the inertial force. (This can be contrasted with the 
cylindrical solution a, where the balance is provided by the cen- 
trifugal force: f B ± > fE±+fB± ~ -fc_L > -f/J_; see Figure 4b). 
As in the nonrelativistic, cold MHD flow considered by Bland- 
ford & Payne (1982), the curvature radius is nonnegligible and 
the solution terminates at a finite height above the disk. 

5. DISCUSSION AND CONCLUSION 

We regard the fiducial solution a as providing the most plau- 
sible representation of an accelerating GRB outflow. Neutrino 
energy deposition and magnetic field dissipation are likely to 
give rise to a nonnegligible thermal component at the origin, 
as envisioned in the original fireball scenario. Such a compo- 
nent is incorporated into solution a but is missing in solution b. 
It is, nevertheless, worth reemphasizing that the dominant en- 
ergy source even in solution a is the Poynting flux. We note 
in this connection that one of the potential problems of tra- 
ditional fireball models has been the "baryon contamination" 

16 As we argue in Vlahakis, Peng, & Konigl (2003), the baryon loading issue may 



issue: Why does the highly super-Eddington luminosity that 
drives the fireball not lead to a more significant mass loading of 
the flow and thereby keep comparatively low? If most of the 
requisite energy were injected in a nonradiative form, then this 
problem would be alleviated. 16 The expectation that GRB out- 
flows are Poynting flux-dominated also renders the slow-rotator 
(enthalpy flux-dominated) solution c less relevant to their inter- 
pretation than the fiducial solution. 

GRB outflows are inferred to be highly collimated at large 
distances from the origin: this has motivated us to construct our 
fiducial solution in the current-carrying (J\\ < 0) regime, which 
is applicable near the symmetry axis (see Fig. 1). However, as 
demonstrated in Figure 3, the illustrative current-carrying solu- 
tion d also becomes well collimated on a similar length scale, 
so it may also provide an adequate representation of such out- 
flows. In fact, this regime may be a natural choice for modeling 
the far-downstream region of the flow since it corresponds to 
a plausible current-configuration on large scales; furthermore, 
as shown in Vlahakis (2003), it also results in a more efficient 
acceleration. We have opted to focus on the current-carrying 
regime in this paper since the solutions in this case, unlike the 
return-current solutions formally extends to infinity. We incor- 
porate the return-current regime into a GRB outflow model in 
Vlahakis at al. (2003). 

Our fiducial self-similar "hot" solution as well as our "cold" 
relativistic MHD solution (and the corresponding ones derived 
by Li et al. 1992 and Contopoulos 1994) have cylindrical 
asymptotics. In the nonrelativistic theory, it was shown (Vla- 
hakis & Tsinganos 1998) that it is possible to derive exact so- 
lutions of radially self-similar flows that have either one of the 
following three types of asymptotes: cylindrical, parabolic, or 
conical. Although we expect that it should be possible to ob- 
tain different asymptotic behaviors also in radially self-similar 
relativistic flows, this has not yet been studied in detail. Our ba- 
sic conclusions about the magnetic acceleration of the flow do 
not, however, depend on the shape of the asymptotic fiowlines, 
although the quantitative constraints on the mass loading (to as- 
sure that the flow is optically thin in the shell-collision region) 
would be eased if the flow continued to expand rather than col- 
limated to a cylinder. GRB jets are often modeled in terms of 
conical jets, and this picture has gained support from observa- 
tions of a panchromatic break in the light curve of several GRB 
afterglows (from which jet half-opening angles in the range of 
2° to 20° have been inferred; e.g., Panaitescu & Kumar 2002). 
However, a sharp physical boundary may not be required for 
interpreting the data in the context of a beam-patterned out- 
flow (e.g., Rossi, Lazzatti, & Rees 2002; Lyutikov & Blandford 
2002). Furthermore, it has been argued that a cylindrical jet 
model could provide a better fit to the light curves of at least 
some afterglows and might perhaps even be able to account 
for the panchromatic breaks commonly attributed to a finite jet 
opening angle (e.g., Dar 1998; Cheng, Huang, & Lu 2001). 

Even if GRB jets are not characterized by a global cylindri- 
cal geometry, the asymptotically cylindrical solutions that we 
derived may be applicable to the "patchy shell" scenario for 
variable GRB outflows (Kumar & Piran 2000). In this picture, 
the outflowing shells are ejected within a given opening angle 
but do not fill the entire solid angle into which they are injected. 
In the context of the internal shock model, this scenario allevi- 
ates the energy requirements on bright bursts; it is also con- 
sistent with the apparent lack of a strong correlation between 

in principle be completely resolved in the context of the MHD acceleration model. 
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the prompt high-energy and the afterglow fluences, and it may 
account for the early-time variability in the afterglow emission 
of a source like GRB 021004 (Nakar, Piran, & Granot 2003). 
In the phenomenological model considered by Kumar & Pi- 
ran (2000), successive blobs of angular scale 1 /70c are ejected 
randomly within an opening angle of 10°, and the number of 
blobs ejected along a given direction (given that more than one 
is required for the production of internal shocks) is also taken 
to be a random number (distributed uniformly in a predeter- 
mined range). This model was intended to mimic the behavior 
of causally disconnected regions in large-scale ejected shells. 
An alternative physical picture is that of a collection of magnet- 
ically active regions on the surface of the disk that feed the GRB 
outflow. In this picture, the magnetic field associated with each 
"patch" guides successive "mini" shells along roughly congru- 
ent paths, thereby enhancing the efficiency of internal-shock 
production. However, the ejection directions from different 
sites need not be parallel to each other, so some of the beamed 
7-ray emission may miss the observer: this might contribute to 
the appearance of quiescent times between GRB pulses (Nakar 
& Piran 2002a). 

The high (~ 50%) Poynting-to-kinetic energy conversion ef- 
ficiency exhibited by all our fast-rotator solutions (solutions a, 
b, and d) has made it possible for our models to accommo- 
date the internal-shock emission scenario without requiring a 
prohibitively large energy input at the source (see §4). How- 
ever, currently available data on the radiated 7-ray energy E 1 
(e.g., Panaitescu & Kumar 2002) and the outflow kinetic en- 
ergy E K (e.g., Frail et al. 2001), which indicate that they are 
both approximately equal (~ 10 51 ergs), are inconsistent with 
this scenario for any reasonable radiative energy efficiency of 
such shocks (see Nakar & Piran 2002b for further discussion 
of this issue). One possible resolution of this apparent incon- 
sistency is that the prompt high-energy emission arises directly 
from electromagnetic energy dissipation without the interme- 
diate step of Poynting-to-bulk-kinetic energy conversion (e.g., 
Thompson 1994; Smolsky & Usov 2000; Lyutikov & Black- 
man 2001; Drenkhahn & Spruit 2002). Our solutions seem to 
be consistent with this possibility: they imply that about half 
of the injected energy is converted into bulk kinetic energy; if 
the remaining Poynting flux can be efficiently converted into 
high-energy emission, then this would lead to an approximate 
equality between £ 7 and E K . Drenkhahn & Spruit (2002) used 
ideal-MHD equations with parameterized magnetic energy dis- 
sipation to demonstrate (in the context of a strictly radial flow, 
i.e., without taking into account the transfield force-balance 
equation that determines the fieldline shape) that the Poynting- 
to-radiation energy conversion could in principle be as high 
as 50% and that the dissipation may also result in significant 
bulk acceleration. It would be interesting to combine the ap- 
proach taken by these authors (see also Spruit et al. 2001 and 
Drenkhahn 2002) with the one utilized in this paper to consider 
the effects of dissipation in a nonradial flow. 

We now briefly compare the results of our exact solutions 
with some of the previous work on MHD effects in GRB out- 
flows. The incorporation of "disposable" magnetic energy into 
the traditional fireball model was discussed by Meszaros et al. 

(1993) . The behavior that they infer is consistent with our so- 
lutions in the regime where 7 oc 03 (but not beyond it). Usov 

(1994) interpreted the outflows in terms of pulsar winds that 
generate intense electromagnetic waves at the point where the 



ideal-MHD approximation breaks down. He highlighted the 
fact that the Lorentz factors of the particles accelerated by such 
waves could be as high as <~ yu 2 / 3 and thus greatly exceed the 
terminal Lorentz factors (<~ /U 1 / 3 ) attained in the Michel (1969) 
ideal-MHD monopole solution. The acceleration efficiency of 
our collimating (and thus nonmonopolar) ideal-MHD solutions 
is much higher yet: they yield 700 [ijl > yu 2 / 3 . 

Meszaros & Rees (1997) examined a magnetized conical jet 
with a pure electron-positron composition. They deduced (and 
our solutions have confirmed) that the pair-photon fluid is ini- 
tially accelerated by thermal pressure, with the magnetic field 
acting only to guide the flow. They also pointed out that radia- 
tion drag effects allow the pairs and photons to remain coupled 
even beyond the point where the scattering optical depth drops 
below 1 and that at some point pair annihilation freezes out. 
Grimsrud & Wasserman (1998) subsequently showed that the 

surviving pairs carry only a fraction <~ 10 7/ of the initial 
energy, implying that, in the absence of magnetic acceleration, 
this scenario is extremely inefficient. We note, however, that 
if a nonnegligible fraction (1 — 7/&//X) of the total injected en- 
ergy is in magnetic form ([with the ratio of magnetic to pairs- 
plus-photons energy injection rates being (y« _ 7/'£i)/7i£! ~ 7/> 
as assumed by Meszaros & Rees 1997), then, based on our 
"cold" solution (see § 4.2.1), half of the magnetic energy could 
be eventually transformed into kinetic energy of the surviv- 
ing pairs, reaching overall efficiencies ~ (/i-£;;7;)/2/i. 17 For 
7, ~ 1 this yields an efficiency of ~ 25%, which is significantly 
higher than the values estimated by Meszaros & Rees (1997). 
As was, however, discussed by these authors, the resulting out- 
flow would have a huge Lorentz factor and might not by itself 
be able to reproduce the observed properties of real GRBs. 

Although the presentation in this paper and its companion is 
focused on the application to GRBs, our formalism is quite gen- 
eral and the solutions that we present should be relevant also to 
other magnetically driven relativistic outflow sources, notably 
AGNs and microquasars. The possible relationship between the 
different classes of beamed relativistic outflows in astrophysics 
has already been noted before by various authors. One poten- 
tial attraction of adopting a common modeling framework for 
these outflows is that it could shed new light on each of the 
different types of sources. For example, the concept of inter- 
nal shocks was originally proposed in the context of AGN jets 
(Rees 1978), and it has recently been adopted also for mod- 
eling microquasars (Kaiser, Sunyaev, & Spruit 2000). In the 
latter class of objects, there is direct observational evidence (in 
the form of correlated X-ray, infrared, and radio flux variations 
and the appearance of superluminal radio knots) for a causal 
connection between a rapid accretion episode onto the central 
black hole and the ejection of superluminal blobs (e.g., Mirabel 
& Rodriguez 1999). This is the same basic scenario as the one 
commonly adopted (albeit without direct observational support) 
for GRB outflows. A similar correlation between the X-ray 
and radio flux behavior and the appearance of radio-bright su- 
perluminal knots was recently discovered also in an AGN jet 
(Marscher et al. 2002). Although the observed behavior in this 
case does not support the notion that AGNs are a simple scaled- 
up version of microquasars, this might be attributable to a dif- 
ferent origin of the underlying disk instability that induces the 
accretion/ejection episodes (e.g., Menou & Quataert 2001) or 
to different environmental conditions (e.g., Heinz 2002). 

There is now also evidence that the high-energy (GeV-range) 



In this estimate we identify po as the rest-mass density of the final pair population, which represents only a tiny fraction of the initial rest-mass density of pairs. 
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7-ray emission in the blazar class of AGNs originates in super- 
luminal radio knots — which likely represent discrete blobs or 
shocks in a relativistic jet — at a large distance from the origin 
(e.g., Jorstad et al. 2001b). 18 This is essentially the picture en- 
visioned for GRB outflows. The superluminal knots in blazar 
jets such as 3C 345 (e.g., Unwin et al. 1997) and 3C 279 (e.g., 
Wehrle et al. 2001) are inferred to move along distinct curved 
paths, although in the former case it has been argued (Wardle et 
al. 1996) that at least two of the knots follow each other along 
the same trajectory. This behavior is consistent with the "patchy 
shell" model discussed above in connection with GRBs. There 
are also indications that the outflows continue to accelerate on 
large scales: for example, the Lorentz factor of knot C7 in the 
quasar 3C 345 was inferred to increase from <~ 3 to > 10 as it 
moved from a distance of r ~ 3 pc from the galactic nucleus to 
r ~ 20pc (Unwin et al. 1997). A scaled-down version of this 
behavior was found in the radio galaxy NGC 625 1 , where knots 
in the radio jets were inferred to accelerate from <~ 0.13 c at 
r w 0.53 pc to - 0.42c at r w l.Opc (Sudou et al. 2000). Such 
large-scale acceleration is most naturally interpreted in terms 
of a Poynting-dominatedjet model of the type discussed in this 
paper (N. Vlahakis & A. Konigl, in preparation). 

In conclusion, we have derived in this paper the general 
formalism for radially self-similar, relativistic MHD outflows 
and presented exact solutions (obtained by solving the Euler 
and transfield equations and imposing regularity conditions at 
the relevant critical points) to illustrate their basic properties. 
We focused on trans-Alfvenic flows, deferring a discussion of 
super-Alfvenic solutions to the companion paper. We consid- 
ered "hot" and "cold" rapid-rotator flows (the latter reproduc- 
ing previous results by Li et al. 1992 and Contopoulos 1994) as 
well as slow-rotator flows in the current-carrying regime (which 



should apply near the rotation axis), but also described a "hot," 
fast-rotator solution in the return-current regime. In all cases, 
we demonstrated that the Poynting flux injected at the source 
can be transformed with high (~ 50%) efficiency into kinetic 
energy of relativistic baryons at a large (but finite) distance from 
the origin. We concentrated on applications to GRB outflows, 
although we pointed out that similar solutions may describe jet 
acceleration in AGNs and microquasars. In the application to 
GRBs, we presented a particular "hot," fast-rotator outflow as 
a fiducial solution. In this case the outflow is initially accel- 
erated thermally, with the magnetic field acting only to guide 
and collimate it, but subsequently magnetic effects become 
fully dominant and are responsible for the bulk of the accelera- 
tion (which occurs between the classical fast-magnetosonic and 
modified fast-magnetosonic surfaces) and for the final collima- 
tion to cylindrical asymptotics. We stress, however, that this 
solution is only illustrative and is not meant to provide an accu- 
rate description of a typical GRB flow. 19 As we discuss in Paper 
II, super-Alfvenic outflows can also provide plausible represen- 
tations of GRB outflows; their basic properties, however, are 
quite similar to those of trans-Alfvenic flows. We have, fur- 
thermore, concluded that solutions in the return-current regime 
may also be relevant to these flows; we address this possibility 
in greater detail in a separate publication (Vlahakis at al. 2003). 
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12635 and by the U.S. Department of Energy under grant 
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APPENDIX 

A. EQUATIONS IN THE AXISYMMETRIC CASE WITH E^ = 

We may combine equations (13a), (13c), and (13d) to obtain 7 , V^, and as functions of M , x, and G using 

03 x l-M z -x z c 

v = — — B+®n4>= — —6 7=^ — a ( A2 ) 

4tt7Po 4tt7Po 03 G l-M 2 -x 2 ' £l-M 2 -x 2 ' 

6*2 / rP dP \ r P 

4nM 2 Vio Poc 2 J {s , A=const} T-lpoc 2 



Knowing the field-line constants (for each s) 

1/2 



m A (A,s)=\^—j , x A (A,s)=—, fMA,s), <7 M (A,s)=-j^-, q(A,s)= {j^i^I J > ( A4 ) 
we can find the quantities x ,M , G , £ (which in general are functions of A , £, and s) by solving the following system of equations: 

x = x A G, M A = q ^ — — , (A5) 



(€-1) 

the Bernoulli equation 



r-i 



p 2 G 2 (l -M 2 -x 2 ) 2 -x 2 (G 2 -M 2 -x 2 ) 2 _ l + alM^_ /ro^A\ 2 (Ag) 



£ 2 G 2 (l-M 2 -x 2 ) 2 i 2 x 4 

18 Although the bulk Lorentz factors inferred from the apparent superluminal motions in blazars are not as high as those indicated in GRB outflows, they can 
nevertheless be substantial, possibly exceeding 40 in a number of cases (e.g., Jorstad et al. 2001a). 

19 For one thing, the value of in this solution is about an order of magnitude higher than the minimum Lorentz factors typically inferred in real GRBs (e.g., 
Lithwick & Sari 2001). 
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which is obtained after substituting all quantities in the identity 7 2 ^1 - V^/c 2 ^ = l + j 2 V 2 /c 2 using equations (A1,A2) and which : 
the nonrelativistic limit takes the familiar form (after Taylor expanding in 1 / c 2 ) 



v 2 + r p (snB d 



:(M-1)C 2 



2 r-ipo *a 

and the transfield equation (obtained from the component of the momentum equation along - VA) 

x A (A,s) 



din 



x 2 (V 5 A) 2 



®f (A '^-LA(l-M 2 -x 2 ) 
oA 



2± 



(W + 



[i 2 x\A 2 



m\a 2 M 2 G^ \ \-M 2 -x 2 



-M 1 



"2^ 



fi 2 A 2 x 6 A 



'M*"A 



r-i 

■V S A-—V S 
1-G 2 ' 



A 2 x\ 



M 2 a 2 M K 



•V 5 A = 0, 



05- V S A- 
•VA- 



where the operator L = V 2 - (2/ 05)03 - Vis related to the curvature radius 1Z = | VA | (LA - VA • V In | VA /GJ | ) 1 . 
The latter equation can be written as 

f G ± + fr± + tc± + ti± + tp± + (e± + fax = , 



(A7) 



(A8a) 



where the subscript _L denotes the component of a force perpendicular to the poloidal fieldline and pointing toward the axis (along 

n = bx (j> = - VA/ | VA |), and where the individual terms are given by 

fcx = 0, (A8b) 
frx = 0, (A8c) 



f/x=- 



t P ± = 



05 

„2„_tT/2 



sin-0 : 



B 2 , c 2 *^x> 2 1 (G 2 -M 2 -x 2 \ 05 -V, 
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f B ± = 



fig 
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ax A /i 



R 2 



1-G 2 



l-M 2 -x 2 
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V 5 a 



4aG 2 



(A8d) 
(A8e) 

(A8f) 
(A8g) 
(A8h) 



B. EQUATIONS IN THE r SELF-SIMILAR CASE 

In the r self-similar model described in §3, the integrals of motion have the following form: 



OQ\ 2 



OM = <J$ 

with the fieldline constants given by 

Box 2 



where <tq , otQ = const , fj,,XA,q = const , A = 



B C3 2 



I a 2 — Q!q 



F-2 

a 2 



L = c®ox A na 



1/2 



r-l (\itc 2 F 



r-i 



(Bla) 



(Bib) 



Note that <jm is the "magnetization parameter" in the monopole solution of Michel (1969). 

In cases where F > 0, the requirement that the magnetic flux function A vanishes on the rotation axis a = implies that chq = 

when um = const and A = ^^a? . We examine solutions with ao = throughout this paper. 
Equations (A5) and (A6) take the form 

x = x A G, (B2a) 



M 



■q 



(£-1)^ ' 



(B2b) 
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// 



2 G 4 (l-M 2 



-x 2 G 2 -M 



2 -x 2 ) 



G 4 (l-M 2 -x 2 Y 



= 1 + 



F 2 a 2 M M A sin 2 f 



£ 2 x 4 cos 2 W> + 60 ' 



and can be solved for x , £, and ip. 

Using tantp = (dz/dXB) A and equation (17), we get tan?/; = d(G/ tan9)/dG, or, 

dG 2 2G 2 cosV' 



d9 sin 9 cos (tp + 6) 



The transfield equation (A7) becomes 



tan(V> + #) 



\-M 2 -x 2 



(F-l) 



4 2 2 



1-G 2 



l-M 2 -x 2 



-sin 2 9 



M 2 + Fx 2 -F+\ 



4 2 2 
X A \X X 

~F 2 a 2 M M 2 



G 2_ M 2_ X , 

\-M 2 -x 2 



+ 2 



cos 2 (ip+9) 

r-i f-2 £(£-i)x 4 



r f 2 c 



M 



M 2 



The magnetization function o is defined as the Poynting-to-matter energy flux ratio: 

-C3fiB /* A c 2 



2 2 

jd -x 



G = 



l-M 2 -x 2 ' 



At the Alfven singular point 8 = 9a, G = 1, x = x A , M 2 = 1 -x A , £ = £a, ^ = 4>a, g = <ja, and, using l'Hospital's rule, 



a A = 



2x 2 A cos IpA 



p A sin 9 A cos(9a + iPa) ' 



l-M 2 -x 



l-M 2 -x 2 



a A + V 



1-G 2 



l-M 2 -x 2 



a A /x 



X G 2 -M 2 



a A +l ' 



l-M 2 -x 2 



x A-( l - X A) a A 

x 2 (a A +l) 



Substituting the above ratios into equation (B2c) yields ^asa function of tp A , 9a, and a A , 



M 2 = 



(VA+l) 2 



A- [x 2 a-va(1-x 2 a )] _ 
whereas equation (B2e) gives the Alfven regularity condition 

F 2 a 2 M (l-x 2 A )(a A +l) 2 sm9 A 



X A^,A + ' 



F 2 a 2 M (l-x 2 ) 2 sin 2 9 A 



x 2 cos 2 (6U + V>a) 



^ 2 cos 2 (6U + Vm) 
T-l(F-2M A -l)(l-x 2 A ) . 



f-2V "' ^ 2 A1 ^ sing A + 2cos^sin(g A + V>A)^ + ^[(F-l)(l-4)-l] 
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(B2c) 



(B2d) 



(B2e) 



(B3) 



(B4) 



(B5) 



(B6) 



Next we prove the following statement: if one chooses T , F (the model parameters) and specifies seven boundary conditions on a 
cone 8 = const, then it is possible to derive all the other model parameters, and the solution is uniquely determined. 
Proof: Suppose that the quantities Bg, B^„ V r , Vg, V<p, po, and P are given as functions of r along the cone 9 = const; i.e., Bg = -C\r F ~ 2 , 
fi = -C 2 r F ~ 2 , V r /c = C 3 , Vg/c = -C A , V^/c = C 5 , po = C 6 r 2 ( F " 2 >, and P = C ir 2{F ~ 2 \ Then one finds 

7=l/(l-V 2 /c 2 ) 1 / 2 , x = {V 4> -VgB^/Bg)/c, </; = ^-0-arctan(-V r /V ), £ = 1 + J\ f , , (B7a) 



r-ipoc 2 ' 

M 2 = 47rp Q aiVg/Bg) 2 , < ? = M 2 (e-l)r L r/C, &m = -~f£,x 2 Vg/ cFM 2 sin 9 , p = ^- f^ ,B ° , 

4njpocVg 

^ ^x^f ' G = X/XA ' ^ B ^ F = - r2 ~ FB ^ GF l^ F ~ l ■ 



(B7b) 



(B7c) 



One is thus in a position to start the integration, using equations (B2c), (B2d), and (B2e), and find a solution that uniquely 
corresponds to the seven boundary conditions (or, equivalently, to the parameters CjJ = 1 , • • • ,7) and the model parameters T and 
F. QED. 

Whether the solution actually hits and passes smoothly through any given singular point depends on the choice of the boundary 
conditions. From a physical standpoint, the most robust solutions cross all the three singular points (the modified slow-magnetosonic, 
the Alfven, and the modified fast-magnetosonic points). In this case, the three regularity conditions at the (a priori unknown) positions 
of the singular points impose three relationships among the boundary conditions. 
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C. ALFVEN AND MAGNETOSONIC WAVES IN RELATIVISTIC MHD 



Suppose that we have obtained a solution of the axisymmetric, relativistic, ideal-MHD equations (3)-(7). If we consider localized, 
fast varying, axisymmetric disturbances, then we may assume that the unperturbed solution is uniform and time-independent and 
neglect all its space and time derivatives. We may then look for perturbations having a Fourier dependence exp[/(wf — k-r)] = 
exp[/(w co f co - fc co • r co )], where, by using the Lorentz transformations between the comoving (r co , f co ) and observer's frame (r , t), we 
get w co =j(uj-k- V) ,fc co = k-V [juj/c 2 - (7- 1) k ■ V/V 2 ] 20 

It is more convenient to analyze the disturbances in the (comoving) flow frame. Define a local Cartesian system of coordinates 
(x,y,z) such that B co = B co z, and fc co = k co (zcos9o+xsm6o). After linearizing equations (3)-(7), we may express all the perturbed 
quantities in terms of the perturbation 5V co . After some manipulation, we obtain 



'D n Di 3 
D 2 2 
,013 D 33 



D n = %&in 2 

c z 



where 



D 22 = 




D l3 = -j sin 8 cos 9 , 



cos 8q 



D 33 = 



C l C os 2 6 -^° 



c 2 k 2 n 



where the Alfven speed can be expressed in terms of the corresponding proper speed Ua, which satisfies U\ 



Besides the trivial entropy wave uj 
modes listed below. 



v 2 /(l- v 2 /c 2 ) = 

(which, however, corresponds to u) ^ 0!), the dispersion relation | D |= yields the wave 



• Alfven waves: D22 = 0, or u) co /k co = ±vaCOS0q, corresponding to a displacement SV co normal to the {B co ,k co } plane. 
Transforming the dispersion relation to the observer's frame, we get 
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sponding to a displacement SV co in the { B co , k co } plane. In the observer's frame, we have 

-1 



— Y 

c 2 k 2 J 



u-k-V 



ck 
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uB p 
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= . 



(C2) 



An interesting property of the waves, related to the discussion on critical/singular surfaces of steady-state MHD, is the following: 
If the component of the flow proper velocity along the wavevector equals in magnitude, but is opposite in sign, to the comoving 
proper phase velocity of the wave, then uj = and the wave is static. (The converse is also true.) Thus, 



:0 



& lV --= UC ° /kc ° 

* (l-u 2 Jc 2 k 2 co ) 1/2 



(C3) 



This statement is easily proved by combining w co = 7 (uj - k ■ V) with the Lorentz invariant u 2 - c 2 k 2 Q = u> 2 - c 2 k 2 and solving for oj. 
Equation (C3) is the generalization of the property of nonrelativistic static waves V ■ k + uj co = uj = 0. It is consistent with the result 
that proper speeds are the appropriate generalization of ordinary speeds in relativistic theory (e.g., Konigl 1980). 
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